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ABSTRACT 


This document presents the mathematical specifications of 
Release 4.0 of the Attitude Determination Error Analysis 
System (ADEAS) , which provides a general-purpose linear 
error analysis capability for various spacecraft attitude 
geometries and determination processes. The analytical 
basis of the system is presented, and detailed equations are 
provided for both three-axis-stabilized and spin-stabilized 
attitude sensor models. 
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SECTION 1 - INTRODUCTION 


The Attitude Determination Error Analysis System (ADEAS) 
provides a general-purpose linear error analysis capability 
for various spacecraft attitude determination processes. 
ADEAS does not process sensor data but simulates the atti- 
tude determination logic and computes the resulting attitude 
determination accuracy. The spacecraft attitude determina- 
tion scenarios that can be analyzed by ADEAS are described 
below: 

0 From low— altitude Earth orbits to International 

Sun-Earth Explorer (ISEE)-3 type of Earth-Sun li- 
bration point orbits 

• Spin— st abi lized or three— axis— stabilized spacecraft 
attitudes 

• Batch weighted-least-squares and' sequential filter 
attitude determination methods 

• Sensor complements, which are subsets of Sun sen- 
sors, Earth sensors, star sensors, gyros, and mag- 
netometers . 

These scenarios include most of the existing and anticipated 
Earth satellite attitude determination systems. A possible 
exception is that attitude rate information is assumed 
available for use in the propagation of satellite attitudes 
in the multiframe method. The rate information is usually 
provided by gyros. 

The ADEAS system requirements are presented in Reference 1. 
The detailed mathematical specifications for ADEAS are pre- 
sented in this document. Section 2 presents the mathemati- 
cal formulations of linear error analyses for batch and 
sequential estimators. The formulations are general and not 
limited only to attitude determination systems. Sections 3 
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and 4 describe the attitude and sensor models included in 
the program. Section 5 specifies the reference systems and 
vectors used iii ADEAS. 

Although many of the algorithms have been extensively re- 
vised, this document is based in large part on Reference 2. 


0450 


1-2 



SECTION 2 - FILTER SPECIFICATION 


This section presents the specifications for error analysis 
using batch and sequential filters. While ADEAS deals spe- 
cifically with attitude determination. Section 2.1 presents 
a widely applicable system model, and Sections 2.2, 2.3, and 
2.4 derive general error analysis equations for batch and 
Kalman filters. Sections 3 and 4 then deal with specializing 
the general specifications to the cases of spin-stabilized 
and three-axis-stabilized spacecraft. 

2 . 1 SYSTEM MODEL 

Let x be an N-dimensional vector that characterizes the sys- 
tem under consideration. This state vector evolves in time 
according to the following dynamic model: 


x(t) = f(x(t), t) + u(t) (2.1-1) 

where the dynamic noise u(t) is a Gaussian white noise proc- 
ess with mean and covariance given by 


E [u( t) ] = 0 

E[u(t) u^t*)] = Q u 6 ( t - f) 

(E [ . . . ] denotes taking the expectation value.) x includes 
all parameters of interest necessary to compute x even though 
some parameters may have zero derivative. For spacecraft 
attitude determination, x includes the spacecraft attitude 
parameters and additional dynamic parameters such as gyro- 
scope biases and alignments. 
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The true value of the state vector is never exactly known 
but is instead estimated by the state estimate vector, x* . 
This estimate evolves in time according to 


x*(t) = f ( x* ( t ) , t) 
The state error vector given by 


( 2 . 1 - 2 ) 


Ax( t) = x ( t ) - x*(t) 

is assumed to always remain small, so linear error analysis 
techniques may be u^ed. To first order, then 


Ax( t ) = x(t) - x*(t) = f(x(t), t) - f(x*(t), t) + u(t) 

» 

(2.1-3) 

(t) Ax( t) + u(t) 

9x 

Integrating this formally gives 


Ax( t) = <J>(t, t Q ) Ax(t Q ) + ( t , t Q ) (2.1-4) 

where the state transition matrix <|> is given by 

i(t, t n ) = *£(t) 4>( t , t ) 

0 3x 0 (2.1-5) 

♦<V t Q ) = I 
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and the random excitation vector i|» by 


Y(t, t Q ) 



It follows from Equation 
property 


4><t 2 , t Q ) 


4>(t, f) u ( t ' ) dt* (2.1-6) 


(2.1-5) that $ obeys the group 
= <i>(t 2 , t x ) Wt^ t Q ) (2.1-7) 


The random excitation vector satisfies the equation 


y(t, t Q ) = <t>(t, t ± ) v(t lt t Q ) + v(t, t x ) (2.1-8) 

A filter produces state estimates based oh information ob- 
tained from measurements made at discrete times. Let y^ be 
a measurement value obtained at time t^. In ADEAS, measure- 
ments performed simultaneously are treated as independent 
scalar measurements, so the times t^ need not be distinct. 
Measurements are related to the state vector by the following 
measurement model: . 


= g i (x(t i ) , p) + v i 


(2.1-9) 
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where p is a vector of measurement parameters and is a _ 

Gaussian white noise process with mean and covariance given 
by 


E[v t ] = 0 



E v. v. =0 for i ^ j 
i 1 


( 2 . 1 - 10 ) 


Note that we will consider p to include the parameters neces- 
sary for all possible measurements, not just those measure- 
ments made at any specific time t^. For spacecraft attitude 
determination, p would include all sensor alignments, biases, 
scale factors, etc. 


The functions g^ are assumed to be known functions of impre- 
cisely known arguments. Therefore, it is possible to com- 
pute expected measurement values by 


y* = g i (x*(t i ), P*) (2.1-11) 

where p* is a vector of estimated measurement parameters. 

The measurement residual between the actual and computed 
measurements is then 


Ay^ = y^^ - y? = g i (x(t i ), p) - g i (x*(t i ), p*) + 

( 2 . 1 - 12 ) 

3g- _ 3g. _ 

~~zr~ Ax ( t . ) + - — Ap + v . 

3x 1 3p 1 
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where Ax and Ap = p - p* are both assumed small. 

The function g i (x(t i ) / p) can be written as 

P) - (1 + kj) hjUCtj), i') + b i (2 . 1 _ 13) 

+ a i sin(to i + t^) 

where k. = measurement scale factor error 
b^ = sensor bias 

ai = amplitude of orbit-related or other unknown 
periodic error source 

o)^ = frequency of periodic error source 

= phase angle of periodic error source 

The parameters k. , b. , a^, and are assumed to depend 

only on the measurement type and not on the measurement time. 
The frequencies and phases i|f^ are assumed to be exactly 
known . 

The parameters k. , b^ and a^ are a subset of the vector of 
measurement parameters; the vector p' contains the remaining 
measurement parameters. The partial derivatives of the 
measurement with respect to these parameters are especially 
simple: 


, hjUCt.), p-> ~ 



( 2 . 1 - 14 ) 


( 2 . 1 - 15 ) 
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(2.1-16) 



sin (to^t^ + 


The partial derivatives with respect to x and p' are derived 
in Sections 3.3 and 4.3. In these sections, they are denoted 
by 8y*/3x and 3y*/3p rather than 3g/3x and 3g/3p. 

2.2 FSTTMATTON AND COVARIAN CE ANALYSIS 

It is usually not necessary to estimate all of the state 
parameters. However, it is necessary to provide estimates 
for those measurement parameters that are not exactly known. 
Therefore, a filter should produce estimates for a set of 
solve— for parameters including a subset of the state param- 
eters and a subset of the measurement parameters . The re- 
maining parameters are then consider parameters whose values 
may contain errors that are not reduced during the estima- 
tion process. ‘ 

The state error, measurement parameter error, and random 
excitation vectors and the state transition matrix are thus 
partitioned as follows: 


Ax(t) 


As x (t) 

AG(t) 


(2.2-la) 


Ap 


AS P 

ay 


(2.2-lb) 
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(2.2-lc) 


<Kt, t c ) = 


v ■ 

^s (t ' t o ) 


*=<>>_ 


+ s (t ' 

t 0 ) j e<t, 
1 

■ V 

0 

: V 1 - V 


(2.2-ld) 


where As(t) = 


4s x (t) 


As 


P J 


= solve-for parameter vector 


A(3(t) = dynamic consider parameter error vector 

Ay = measurement consider parameter error vector 

The error equations (2.1-4) and (2 . 1-12 can then be re- 
written as 


As ( t) = <X>(t, t Q ) As(t Q ) + 0(t, t Q ) AfJ(t Q ) + T(t, t Q ) (2.2-2) 


A(3(t) = <t> B (t, t Q ) AR(t 0 ) + »|/ 3 (t, t Q ) 


(2.2-3) 


Ay i = G i As(t i ) + T i Ay + 


(2.2-4) 


where 


4>(t, t Q ) = 


4> s (t, t Q ) ! O 


(2.2-5) 
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( 2 . 2 - 6 ) 


*(t, t 0 > = 


♦s (t ' t o ) 



e(t. 

0 

3g. 

3g ^ 

9s 

_ 9i x 

r = 

3g. 


3g t 

3s_ 


3y 


(2.2-7) 


( 2 . 2 - 8 ) 


(2.2-9) 


!£i 

3x 


ag'i 

i 

i 

! o 

8sL 

1 

1 

X 

1 


( 2 . 2 - 10 ) 


3^ 

3p 


3 g ^ i 3 g ^ 


0 Sp j 3 y 


( 2 . 2 - 11 ) 


We have assumed that no measurements depend directly on any 
dynamic parameter that is not a solve— for parameter. Fur- 
ther, we have assumed that the time evolution of the dynamic 
consider parameters does not depend on any of the dynamic 
solve-for parameters. 

The function of a full estimation system is to determine an 
estimate s* (t) given measurements y i . ADEAS, however, does 
not actually compute an estimate but determines how good an 
estimate would be if it were produced in a given situation. 
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ADEAS does this by computing the estimation covariance ma- 
trix defined by 


P(t) 5 E [ As ( t ) As T (t)] (2.2-12) 

ADEAS thus performs linear covariance analysis for batch and 
sequential filters. The covariance matrix P(t) then pro- 
vides a statistical measure of how good an estimate could be 
produced at time t of a given scenario. 

The random excitation enters into this computation in the 
form of the random excitation covariance matrix, which is 
defined as 


d(t, t Q ) 


= E t 0 > t* (t< t 0 )j 


a s (t ' ‘o' 


a sli^‘' ‘o' 


a sl5<‘' ‘o' 


a B (t, t 0 ) 


where 


(2.2-13) 


a s<‘- ‘o' 


= e [y s (‘' ‘o' ‘ 0 >] 

t 

= f [<J> s (t, t*) i e(t, t-)] Q u 


4>g(t, f) 


e (t, t’) 


(2.2-14) 
dt 1 


d B (t, t Q ) 


= E [V fc ' t 0 ) *B (t ' t 0 > ] 

--J [° i V‘- t- >] °u 


|4> B (t, f) 


(2.2-15) 
dt • 
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(2.2-16) 


d s(3 (t ' t 0 ) = E [*8<t. t 0 ) ^3 (t ' t 0>. 


A . 

1 

1 

o 

II 

* <t, f> ! 9(t, t') 

L 


4>J(t, f) 


2.3 BATCH FILTER 

A batch filter produces an estimate s*(t Q ) at an epoch time 
t , based on a single batch of measurements y that may have 
been made at various times. Thus, 


y 




(2.3-1) 


where each y^^ is a scalar measurement. Similarly, 


and 



(2.3-2) 


The batch filter produces an estimate s*(t Q ) that gives y*, 
which minimizes the cost function 


V 


Ay 1 WAy + As* T W A As* 


(2.3-3) 
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where W E positive-definite symmetric measurement weight 
matrix 

As* = s * - s*(t Q ) = s* - s7t 0 ) + s7t o ) - s*(t Q ) 

= Al(t c ) - As A 

s* E a priori estimate of s*(t Q ) 

W E non-negative-definite synunetric a priori weight 
A matrix 

2.3.1 ESTIMATION ERRORS 

Since the batch filter determines s*(t Q ), it is necessary to 
relate Ay to As(t^) . Substituting Eguation (2.2—2) into 
Equation (2.2-4) gives 


Ay t = G.[0(t i , t Q ) As(t 0 ) + ©(t^ t Q ) Aft(t o ) + t Q )] 

+ r. + v. (2.3-4) 

= F i Al(t o ) + B. A(3(t 0 ) + T i Ly i + 
where, using Equations (2.2-5) through (2.2-8) 



•r 

O' 

CO 

I 


F i = <3i *<*=!• t Q ) = 

— M (t. , t ) i 

— S V 1 O • 

L 3s x 

3i p. 


B . E 

l 


G. 


l 


0(t i , 


t o ) 



©(t^ 


t ) 
o 


(2.3-6) 
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(2.3-7) 


Then, 


where 


u i 5 G i '■'‘V V * W t o ) 


Ay = FAs(t Q ) + Ae 


Ae E BAB ( t ) + TAy + U + N 



” F l" 


"V 


_r r 

F = 

• 

• 

B = 

• 

• 

r = 

• 


• 

“f~ 


• 

“b” 


• 

~r~ 


m_ 


mj 


m 


u = [ U 1 U m] 

T 

N = [v 1# . . . , vj 


(2.3-8) 


(2.3-9) 
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Substituting Equation (2.3-8) into Equation (2.3-3) for the 
loss function then gives 


V = [FAs(t Q ) + Ae] T W[FAs(t 0 ) + Ae] 


+ As(t ) - As, 

O i 


W A Us(t 0 ) " AS A ] 


= As T (t Q ) |w A + F T Wf| As ( t Q ) + As T (t o )(F T WAe - W^As A j 
+ (a^WF - As^W A j As ( t Q ) + A^WAe + As^W a As a 


j^s<t 0 ) + W-^F 1 WAe - w A is A )j T W N [As(t 0 ) 

+ wij 1 (p T wie - w a a; a )J 
- (fW - w a aS a ) k-^f 1 wAi - w a aI a ) 


(2.3-10) 


+ Ae^Ae + As^W a As a 


where 


W„, = W, + F T WF = normal matrix (2.3-11) 

N A 

The final equality in Equation (2.3-10) is valid as long as 
W N is nonsingular. The singularity (or ill-conditioning) 
of W„ indicates a lack of observability of the solve-for 
parameters from the measurements y. 
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If W is nonsingular and positive-definite, then it is clear 
from the form of Equation (2.3-10) that V is minimized when 


As(t c ) = -W' 1 (F T W4e - w A is A ) 

- „ - ) (2.3-12) 

) + TAy + U + N] - W a As' 


= -W N 1 jF T W[BA(3(t 0 , 


= ^S n (t 0 ) + As y (t 0 ) + As R (t 0 ) + As u (t 0 ) 


where 


As n (t 0 ) = W" 1 (W a As a - F T WN) (2.3-13) 

As (t ) = -W" 1 F T wrAY (2.3-14) 

Y o N 

Ai [3 (t 0 ) = -W^ 1 F T WBAB(t Q ) (2.3-15) 


As (t ) = -W., 1 F T WU (2.3-16) 

U O N 

The batch filter produces an estimate s*(t Q ) at the epoch 
time t . This estimate may then be propagated to any other 
time t using Equation (2.1-2). In doing this the estimation 
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errors As propagate according to Equation (2.2-2). Substi- 
tuting Equation (2.3-12) for As(t Q ) into Equation (2.2-2) 
gives 


As ( t) = <fr(t, t c )^As n (t 0 ) + As y (t 0 ) + As 0 (t o ) + As u (t 0 )J 

+ ©(t, t Q ) AB(t 0 ) + Y(t, t Q ) (2.3-17) 

= As n (t) + As^,(t) + AS|j( t) + As u (t) 

where 

As n (t) = 4>(t, t Q ) As n (t Q ) (2.3-18) 

As Y (t) = <t>(t, t Q ) As Y (t D ) (2.3-19) 

ASg(t) = 4>(t, t Q ) Asp (t Q ) + 0(t, t Q ) AB(t 0 ) (2.3-20) 

As u (t) = <D(t, t Q ) Ai u (t 0 ) + Y(t, t Q ) (2.3-21) 

2.3.2 COVARIANCE 

Equation (2.3-17) gives the estimation errors induced by the 
systematic error sources Ay and AR, the random error sources 
U and N and the a priori error As"*. It is assumed that all 
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these error sources 
P(t) is 


are uncorrelated. The covariance matrix 


P(t) = E[As(t) As^t)] (2.3-22) 

= P n (t) + P y (t> + P 0 (t) + P u (t) 

where 

P n (t) = E[As n (t) As^(t)] (2.3-23) 

P^t) E E [ As^t t ) As^(t)] (2.3-24) 

P 3 (t) = E[As B (t) As^(t)V (2.3-25) 

P u (t) = E[As u (t) As^(t)] (2.3-26) 

The following subsections discuss, in turn, each of these 
contributions to the overall covariance. 

2.3.2. 1 Data Noi se Contribution 

From Equations (2.3-18) and (2.3-23) we have 

P n (t) = <l>(t, t Q ) P n (t Q ) * T (t, t Q ) (2.3-27) 

where from Equation (2.3-13) 

P n (t O ) = W N 1 ( W A P A W A + pTwRWF ) W N X (2.3-28) 
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with 




P A = E [ AS A Ai l] 


R = E [NN 7 ] = diag 

\ 

The a priori weight matrix is given by 


W. = p, 1 

A A 


where the a priori covariance is 


P A = 


P Aattitude j 
* 

0 

0 1 

P Aother _ 


(2.3 


(2.3- 


(2.3- 


with P being the a priori covariance of the att 

tude error parameters, as given in Section 3.2 or 4.2.1, 

p being the a priori covariance of other solve-for 

Aother 

parameters, assumed to be diagonal: 


P Aother = diag (a!?,, a? ) (2 .3 

AJ - other 


W is assumed to be a diagonal matrix of the form 


W = diag 


( 


w 


1 ' 


w 


m 


( 2 . 
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30a) 


30b) 


and 


30c) 
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Equations (2.3-28) through (2.3-31) give 


W 



w. 


+ 


E 



F. 

l 




(2.3-32) 


where 


m 



i = l 


(2.3-33) 


with the row vector F. given by Equation (2.3-5). For mini- 

1 _i 

mum variance weighting, set W = R so that P (t ) assumes 

1 n o 

its minimum value P (t ) = W~ . Note that this only 

n o N 

minimizes the data noise contribution to^the total covari- 
ance . 

2. 3. 2. 2 Consider Parameter Contribution 

It is assumed that all consider parameters are uncorrelated 
so that 


E[AyAy ] = diag 



E 


A0(t o ) A[B T (t o ) 


diag 




(2.3-34) 
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Using this with Equations (2.3-14), (2.3-15), (2.3 19) 

(2.3-20), (2.3-24) and (2.3-25) gives 


n 


P y (t) 


= ♦<*' t o ) g 


&S-(t ) 

3Y^o ; 


(t ) 

o> 


<t> T (t, t c 

( 2.2 


n 


3 


p B< t> - Z a B. 

i=l 


t D ) ^-<t 0 ) +e i (t. t Q > 

33^ 


J (2 .: 


‘o' ftr^o 5 + e i (t - ‘o' 


where, from Equations (2.3-14) and (2.3-15), 


m 


* 

3y 


-"n 1 Z F i r l“i = 


i = l 


5jS_/ t ) S 

_3y 1 ( o ) | 


* j 9y n <V 
y 


( 2 .: 


m 


3s 


af<‘o' - -V Z F I B i M i = 


i = l 


( t ) 

L 8B 1 ° 


’ * j 33 ^ C o ; 
’ n 3 


( 2 .: 


with 


0(t, t Q ) = 


0 x (t, t 0 ) ! 


i ! 


! % (t ' V 


) 

-35) 


-36) 


-37) 


-28) 
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I\, and B i defined by Equations (2.2-9), (2.3-5), and 

(2.3-6), respectively; and and W N as in Sec- 
tion 2. 3. 2.1. The quantities and are row vectors of 
dimensionality n^ and n^, respectively. 

2. 3. 2. 3 Dynamic Noise Contribution 

From Equations (2.3-16), (2.3-21), and (2.3-26) we have 


P u (t) - D(t, t 0 > + 4>(t, t 0 ) P u (t c ) 4> T <t, t 0 ) 

(2.3-39) 

- <t>(t, t Q ) w~ 1 f t we [u ^(t, t Q )j 

- Ep(t, t Q ) U 1 ] WFW" 1 <t> T (t, t Q ) 


where 

D(t, t Q ) = E^(t, t Q ) ¥ (t, t Q ) 

with d defined by Equation (2.2-14), and 
s 

p u (t 0 ) = W“ 1 F T WE[UU T ] WFW- 1 (2.3-41) 

Using Equations (2.3-7), (2.3-9), and (2.3-31) gives 


d s (t ' t o > j- 

4 — 


! 0 


(2.3-40) 


m 


-E 


F T WU . > M f (t i , t ) 
i-1 1 


( 2 . 3 - 42 ) 
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where 



(2.3-43) 


with defined by Equation (2.3-5), and 9gV3s x defined by 
Equation (2.2-8). 

The row vectors 9g./9s have dimensionality equal to the 
number of dynamic solve-for parameters. 

From Equation (2.3-42) 


FWUUWF = 


m 


X M G. W t o> 

i = l 1 


m 


X ”g- W t ° ) 


j=l 


m 


- M G. W t O ) *S (t i' t O > M G i 

* . J * 


i=l 


+ M 


(2.3-44) 


G. 

l 


3-1 


X W t o> *S<V V M ' 


j=l 


G. 

3 


i-1 


X m G. W t o > V 


j=l 


M 


G i, 


2-21 


0450 



Using this in Equation (2.3-41) 


m 


w 


- V £ \% w v N l t 

i = l ( 1 


M 


G . 


i-1 


S a s (t i' V 


W, 


l Li-i 

i-l 


(2.3-45) 


^2 M G. d s T(t i / t j ) 


L'j =1 


m t I w d 
n 


where d g is defined by Equation (2.2-14) and for t i > tj 

d’Ctj, t.) = s[? s (t i( t 0 ) ^(t^.t 0 )] (2.3-46) 

Using Equation (2.2-lc), we see that this is the upper left- 
hand corner of the partitioned form of the larger matrix 


d ' ( t i , tj) 


E [‘*'(t i , t Q ) ^(t y t Q )] 

a s< V v I a sB (t i- V 
d'£<V t ’ > * a » (t -' t ' ) 


I3 V i 


(2.3-47) 


This can be written, using Equation (2.1-8) and the fact that 

E [^(t i , tj) «|» T (tj , t Q )] = 0 


for t^ > t j , as 


d ’ ( t^ , tj) = 4> ( t i , tj) d( t j , t Q ) 


l -1. 


(2.3-48) 


= *(t., t Q ) <|> x (t jf t Q ) d(t j , t Q ) 
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The group property of <t>, Equation (2.1-7) has been used to obtain 
the last equality. Then combining Equations (2.2-ld), 

(2.2-13), (2.3-47), and (2.3-48) gives 


a-CV tj ) 


-i, 


= t<t> s <t i . t 0 ) | ettj, t 0 >] 4> <t j( t o > 


(2.3-49) 


Substituting Equation (2.3-49) into Equation (2.3 45) then 
gives 


m 


i = l 


W - "n 1 } S "GjW V n l i + "<3^1 + Q i M Gi 


w 1 

N 

(2.3-50) 


where 


i-1 

Qi = £ *'hty t Q ) 

j=l 


d s (t j' fc o>. 


L d sG (t j' t o ) 


Mi, (2.3-51) 

G j 


and 


Ml 5 M 
G i G i 


♦s (t i' t o ) I 0(t i' t o>‘ 


f t 

1 


a g 

7^ W V 

OS 


B. 

l 


(2.3-52) 


w. 


where is defined by Equation (2.3-6), and Equa- 
tion (2.3-43)- has been used in the last step. 
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From Equations (2.3—42), (2.3-46), (2.3-49), and (2.3-51) we 
also have: 


T 

F WE 


[u ^(t, t 0 >] - £ M^E [-|. 8 (t i . t c ) ^(t, t 0 )] 


i = l 
k-1 


- X) M Gi E [ W V *l lt ‘ t o\ 


m 


S "g.WV V t o ) 

i=k 1 L 


(2.3-53) 


k-1 


m 


- X) m g. d s T(t - V ♦ X) M G. a s<V fc > 

1 i=k 1 


i = l 


= Q 


v' 

+ Q,* ♦ 1 (t, t ) 

V*' V 

1 

CD 

H 

rr 

% 

n 

o 

l 

k x o 

a Il3 (t ' V 


where k is chosen so that t^_ 1 < t < t k and 


m 


k " S M G i 

i=k 


(2.3-54) 


Then, 


T 

F WE 


u r^t, t Q >] = r f t we | u ^(t, t 0 )j 


(2.3-55) 


and P u (t) can be computed by substituting Equations (2.3-40), 
(2.3-45), and (2.3-55) into Equation (2.3-39). 
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2.3.3 COMPUTAT ION 


The analysis equations presented in the previous subsections 
have been carefully derived to allow efficient evaluation. 
This subsection outlines a procedure to efficiently perform 
these computations. This procedure is only intended to show 
the overall structure of the computations and does not neces- 
sarily cover all microefficiency details. The procedure has 
two stages: it first computes analysis results at the epoch 

time t , and then uses these results to compute the various 
covariance matrices at a set of "output times" 

It is assumed that t Q < t.^ < t 2 < ... < t m and that t Q = 

< t 2 * * * < T 2,’ 

Stage 1: Compute Epoch Errors 


The stage 1 procedure assumes that the following are avail- 
able: 


a . 

l 

See 

Equation 

(2.3-29) 

W A' W i 

See 

Equations 

(2.3-3) 

3g . /3s 

*1 X 

See 

Equation 

(2.2-8) 


4.3 



F. 

l 

See 

Equation 

(2.3-5) 

F i 

See 

Equation 

(2.2-9) 

4.3 



B. 

l 

See 

Equation 

(2.3-6) 

D(ti, t 0 ) 

See 

Equation 

(2.3-40) 


4.2 



V 

See 

Equation 

(2.2-5) 


4.2 



W t o > 

See 

4.2 

Equation 

(2.2-ld) 
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9(t., t ) See Equation (2.2-ld) and Sections 3.2 and 
1 o 

4.2 

4> -1 (t i , t Q ) See Equations (2.1-5) and (2.1-7) and Sec- 
tions 3.2 and 4.2 

The procedure produces the covariance matrices P R and P u and 
the sensitivity matrices 9 s/9y and 9s/9B, all at the epoch 
time t Q , and the intermediate matrices , P N , Q', and 

Q . These results are accumulated during a single pass over 
F 

all the measurements. 


1. 

Initialize M , M , M, 
n u 

2. 

For 

i «- 1 to m, do 


a . 

M „ - M n * F I F i 


b. 

- m y . f* r 4 

T ^ g i 


c. 

m g * F l a - w i 


d. 

m f *- F I F i "i 


e . 

M b - F* Bi w. 


f . 

Let M f - |-M Fg j 
sider columns) 


g- 

M u - M u + M g d s 


h. 

M fl - M f) + M B 


i . 

Q. , 4- Q. + <J> -1 
v i+l v l 


j • 

Q f * q f + m f 


k. 

Q* Q* + M g 

i 


2 2 


[ "I (with solve-for and con- 

F pJ r 1 

«b. - [ m f s i m b] 


jtT »» t /-\T TUT ♦ T 


d s (t i' fc o> 


d sB (t i' V 


m; 
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3. 


4. 


N 


n 


u 


<W A ♦ Q p > 


N 


-1 


P N (W a + M„) P 


P„ M P„ 
N u N 


n' 


N 


, for minimum variance 
weighting 

, otherwise 


6 . 


dy 


-P„ M 
N Y 


7. 


J*-’ * " B 


Stage 2; Compute Error Analysis Results 

In addition to the results of stage 1, the stage 2 procedure 
assumes that the following are available: 

a , ff R See Equation (2.3-34) 
y i l 

<J>(ts, t Q ) See Equation (2.2-5) and Sections 3.2 and 4.2 
D ( , t ) See Equation (2.3-40) and Sections 3.2 and 4.2 
9( T i' t o ) See Et 3 uation ( 2 - 2 ~ 7 ) and Sections 3.2 and 4.2 

The procedure produces the covariance matrice P n (T^), 

) , P^x^, p u ( T i) » an{ i P('r i ) for each output time 

1. Set k ■*- 1 


2 . 


For i *• 1 to i, do 


b. 


P n< T i> * *< T i' V P n * < T i' V 


Let ^ = <P(t., t ) %t ) = 
3 y ° 3Y 



- 

3s | 

j 3s 

97 i ! 

' • * j 8 ^ n 

Y _ 


0450 


2-27 



in 


n 


W 


I 

j=i 


2 ls_ lS_ 
Yj SYj 9Yj 


Let ^ = 4 >(t., t ) %t ) + ©(t^ t ) 
9(3 1 ° 3(3 


ls_ 


93. 


_as_ 


93 


n 


3 


m 


n, 


3 


, v ^ 2 9? 

< T i> * Z a 3 • 93, 


3 = 1 


2 <Ls_ 9 s T 
3 -D 83 D 


While t k < x i , do 
(1) Q* <- Q* - M* 


(2) k •*- k + 1 



Let P' = <J>(T i , t Q ) P N j 

[ q t 

'♦S<v t o ) i 0 

e . 

Q k 

m 1 


i 


_e 1 (x i , t o > ! o J 


+ Q ,T 4> _1 (T i/ t Q ) 


d s< T i' V 


d s3< T i' V 


0 

0 


m 


W - D <V V 

+ 4>(Ti, t 0 ) P u 4> T (Ti, to) - P' - P' T 
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P(x i ) «- P n (Ti> + P/V + V T i } + P u (T i ) 


f . 

2.4 SEQUENTIAL FILTER 

A sequential filter produces an estimate s*(t) based on 
measurements taken at discrete times t i < t. Between 
distinct measurement times, the state estimate x*(t) is pro- 
pagated using Equation (2.1-2). For each measurement y i , 
the solve-for parameters are updated based on the pre-update 
state x*(i-) and the measurement. Typically, this update 
has the following form: 


s* ( i+) = s*(i-) + K i Ay i (2.4-1) 

where s*(i+) and s*(i-) denote the estimate of solve-for 
parameters immediately after and immediately before incor- 
porating the information contained in the* measurement . This 
notation must be distinguished from s*(t^), which denotes 
the solve-for parameter estimate incorporating the informa- 
tion contained in all the measurements at t^, which may 
include measurements other than y^. The gain matrix 
determines how much the propagated state is corrected, based 
on the measurement residual Ay^; this is a column vector 
with dimension equal to the number of solve-for parameters. 

2.4.1 ESTIMATION ERRORS 

The estimation error immediately after an update is 


As ( i+) = s(t i ) - s*(i+) = s{t ± ) - s* ( i-) - K i Ay i 


(2.4-2) 


= As(i-) - K^Ay^ , 
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since the true state is continuous at t^. Substituting Equa- 
tion (2.2-4) for Ay^ gives 

As ( i + ) = (I - K i G i ) As(i-) - K^TjAy + v i ) (2.4-3) 

We divide As(t) into separate contributions due to measure- 
ment noise, measurement consider parameters, dynamic con- 
sider parameters and dynamic noise: 


As ( t ) = As n (t) + As^(t) + ASjj(t) + As u (t) (2.4-4) 

These obey the update equations 


A s n (i+) = (I - K i G i ) As n (i-) - «iv i (2.4-5) 

As (i+) = (I - K i G i ) As^( i-) - K-T^y (2.4-6) 

As R (i+) = (I - K.G i ) As (J (i-) (2.4-7) 

As u (i+) = (I - K.G i ) Al u (i-) (2.4-8) 


consistently with Equations (2.4-3) and (2.4-4). If t^=t^ + ^, 
then As ( i+) = As(i+1-), and similar relations hold for As n , 
As^,, Asp, and As u - 

If t. * t i+1 , which means that y^^ is the last measurement 
processed at t^, As(t^) is equal to As(i+), since it 
incorporates the information contained in all the measure- 
ments for times less than or equal to t^. Then from 

2-30 

0450 


Equation (2.2-2), propagating the estimate s*(t._) will 
result in an error for t^ < t < t^ + ^ given by 

As(t) = <I>(t, t i ) As(t i ) + 0(t, t^) + Y(t, t i ) (2.4-9) 

The limit of As(t) as t approaches t i+1 is the estimate 
error As(t. +1 -). Inserting Equations (2.2-3) and (2.4-4) 
into Equation (2.4-9) gives the following propagation equa- 
tions for the noise and consider components of the estimation 

error 


As n (t) = <t>(t, t t ) As n (t i ) (2.4-10) 

As y (t) = <I>(t, t ± ) As^(trr) (2.4-11) 

is R (t) = $(t, t.) As | .(t i ) + 0(t, t^ ^(t^ t Q ) AR(t Q ) 

15 (2.4-12) 


As (t) = ^(t, t.) As u (t^) + ©(t, t^) ^(t^/ t 0 > 

(2.4-13) 


+ ¥(t, t i ) 


The complete specification of As n , As^, As^, and As u requires 
the initial conditions 


As n ( V * is(t o> 


(2.4-14) 


AS (t Q ) > As B (t 0 ) = As u (t 0 ) = 0 


(2.4-15) 
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2.4.2 COVARIANCE 


Equation (2.4-4) gives the estimation errors induced by the 
systematic error sources Ay and AG and the random error 
sources u and It is* assumed that all these error 

sources are uncorrelated. The covariance matrix P(t) is then 


P(t) = E [ As ( t ) As^Ct)] 

= P„(t) + P y (t) + P R <t) + P u (t) 


(2.4-16) 


where 


= 


As n (t) 


As^(t) 


+ D*(t, 


fc o> 


(2.4-17) 


P y (t) = E|Al Y (t) As^(t)J (2.4-18) 

P R (t) = E^As B (t) AsjJ(t)j (2.4-19) 

P u (t) = E[Ai u (t) Al^(t)] - D* ( t / t Q ) (2.4-20) 

The matrix D*(t, t ) represents an estimate of the dynamic 
noise used to compute the gain matrix. It obeys the update 
equation at measurement y. : 

D*(i+, t Q ) = (I - K i G i ) D* ( i- / t Q ) (I - K i G i ) T (2.4-21) 
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If * t i+1 / D* ( t / t Q > obeys the propagation equation for 
t. < t < t. 


•i+1 


D*(t, t Q ) = <t»(t , t i ) D*^, t Q ) <I> (t, t i ) 


( 2 . 4-22a) 


+ D*(t, t i ) 


with 


D*(t, t.) 



(2 . 4-22b) 


and 


d|(t, tj) 


■/ 


[<t> s (t, f) 


e(t, t*)] q * 


fc i 


4>g(t, f) 


e A (t, f) 


("2 


dt ' 

. 4-22c) 


Q* is a diagonal matrix of estimates of the covariances of 
the dynamic noise. The matrix D*(t, t ) has the initial 
value D* ( t / t Q ) = 0. 

The following subsections discuss each of these contributions 
to the overall covariance. In each case, the propagation 
step is unnecessary if t^ = t^ +1 * 

2. 4. 2.1 Noise-Induced Contribution 

From Equations (2.4-10), (2.4-17), and (2.4-22) we have 

P n (t) = 4*(t , t i ) P n (t i ) <l> T (t, t ± ) + D*(t, t i ) (2.4-23) 
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for t. < t < t^ +1 . From Equations (2.4-5), (2.4-17), and 
(2.4-21), the covariance update is 


P n (i+) = (I - K i G i ) P n (i-)(I 


K i G i) T + K i R i K i * (2.4-24) 


with 



(2.4-25) 


The initial value P (t ) is the a priori covariance P A , 
given by Equation (2.3-30b). 

For Kalman filtering, K i is chosen to be the Kalman gain 
defined by 

K i - p n (i -> G i[ G i P n (i -> G i + R i] (2 ’ 4 ' : 

Note that a Kalman filter minimizes only the noise-induced 
contribution to the overall covariance. 

2. 4. 2. 2 Consider Parameter Con tribution 

It is assumed that all consider parameters are (initially) 
uncorrelated so that 


E [ Ay Ay ] = diag 




(2.4-27) 


E 


Af3(t 0 ) AB ( 



diag 




(2.4-28) 
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Using this with Equations (2.4-18) and (2.4-19) and 
8 s/3y and 8s/83 so that 


As (t) = ^E(t) Ay 
Y 8Y 


is B (t) = *5 <t) W3(t > 
15 813 


then gives 


n 


P Y <t) 


i 

-z 

i = l 


8s 

3Yi 


(t) 


8s 

3Yi 


1 T 


(t) 


n 


13 


P 13(t) 


-Z 

i = l 


13. 


Is- 

as. 


(t) 


3s_ 

8S i 


(t) 


where, from Equations (2.4-11) and (2.4-12), 


Is 


Is./ 


±<t) = 4>(t, t.) “(t.) = 
3Y 1 3Y 


M. ' 


! Is 


3ir l j i 3 ^n 


YJ 


8s 


* (t) = 0(t, t.) (t.) + 0(t, t i ) 4> B ( t ± , t Q 

3B 3R 


3s 1 

dti 1 1 


! 8s 


i • • • i 


• i 813 


n 


13 
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defining 


(2.4-29) 


(2.4-30) 


(2.4-31) 


(2.4-32) 


(2.4-33) 


) 

(2.4-34) 



Using Equations (2.4-29) and (2.4-30) in Equations (2.4-6) 
and (2.4-7) gives the following update equations: 


^ (i+) = (I - K.G.) (i-) - K.T. (2.4-35) 

a Y 1 1 3y 


^ (i + ) = (I - K.G.) (i-) (2.4-36) 

8 13 1 1 313 

Both these matrices are zero at the initial time t Q . 

2. 4. 2. 3 Residual Dynamic Noise Contribution 

From Equations (2.4-13), (2.4-20), and (2.4-22) we have 

P u (t) =* f(t, t.) P^t^) o T (t, t L ) + t A ) P u(J (ti) t L ) 

+ e<t. t £ ) pJ (J (t i ) * T (t, t t ) + 0 (t, t L ) t o > e T (t. t t ) (2.4-37) 

+ d( t , t A ) - D*(t, t.) 


where 


p uB (t) = E 


As u (t) 


♦B (t ' t o ) _ 


(2.4-38) 


and dp is defined by Equation (2.2-15). 

It follows from Equations (2.1-8), (2.2-lc), and (2.2-ld) 
that 


V fc ' t o ) “ V fc ' fc i> W t o ) + V*' t i ) (2.4-39) 
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Using this and Equation (2.4-13) in Equation (2.4-38) gives 


P u |3<t> = *<t, t.) P^Ctp +5<t, t.) 

+ ©<t, t t ) a,,^, t 0 > ^(t, tj) (2.4-40) 

+ ° s a^- V 


where 


D sl3 (t - t i > S E [' ! ' <t ' t i > ‘i'] " 


a sB (t ' t i > 


(2.4-41) 


with d „ given by Equation (2.2-16). 
sli 

Using Equation (2.4-39) in Equation (2.2-15) gives 


a„(t, t c ) = + B <t, tp a B (t i , t c ) 4> B (t, tp + a B (t. tp 


(2.4-42) 


From Equations (2.4-8), (2.4—20), and (2.4 21), the 

covariance updates for P u and P u ^ are 

P u (i + > = (I - Kpp P u U-> (I - KjG^ 1 

P u(3< i + ) ■ (I - K i G i ) P uB (i -> 

Both P u and P u ^ are zero at the initial time t Q . 

2.4.3 COMPUTATION 

The analysis equations presented in the previous subsections 
have been carefully derived to allow efficient evaluation. 
This subsection outlines a procedure to efficiently perform 
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these computations. This procedure is only intended to show 
the overall structure of the computation and is not concerned 
with all microefficiency details. The procedure propagates 
an initial covariance P A to the output times x^, . .., x^, 
performing measurement updates at times t ^ , ..., t m . It is 
assumed that t Q < t^ < t 2 < ... < t m and that t Q = x^ < x 2 < 
... < 


procedure 

assumes 

that the 

following are available: 


P A 


See 

Equation 

(2.3-30) 

and 

Sections 3.2 

and 



4.2 






a . 
1 


See 

Equation 

(2.4-25) 




G i 


See 

Equation 

(2.2-8) 

and 

Sections 3.3 

and 



4.3 






r. 

i 


See 

Equation 

(2.2-9) 

and 

Sections 3.3 

and 



4.3 



— * 



o 

rt 

s 

t') 

See 

Equation 

(2.2-5) 

and 

Sections 3.2 

and 



4.2 






e(t. 

t’) 

See 

Equation 

(2.2-7) 

and 

Sections 3.2 

and 



4.2 






4>(t , 

f) 

See 

Equation 

(2.1-5) 

and 

Sections 3.2 

and 



4.2 







f) 

See 

Equation 

(2.2-ld) 

and 

Sections 3.2 

and 


4.2 






D ( t , 

f) 

See 

Equation 

(2.3-40) 

and 

Sections 3.2 

and 



4.2 






d( t , 

f) 

See 

Equation 

(2.2-13) 

and 

Sections 3.2 

and 



4.2 






d B (t, 

f) 

See 

Equation 

(2.2-15) 

and 

Sections 3.2 

and 


4.2 
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D 


sB 


(t. 


f) 


See Equation (2.4-41) and Sections 3.2 and 
4.2 


D* ( t , t') See Equation (2.4-21) 

o , a„ (t ) See Equations (2.4-27) and (2.4-28) 
Y. Bi o' 


The procedure produces the covariance matrices P n (x^), 

P (x.), P rt (x.), P (x.)/ and P(x.) and the analysis matrices 
Y N i 7 13 i u i i 

4>(Xi, T i _ 1 ) and d(x i# x i ;L ) for each output time x i : 

1. Set P„ «- P. 

n A _ 

Initialize 0 s/8y, 8s/ 8B, P u , P u(r d, and d R to zero 
Initialize <}> to the N x N identity matrix 


Set i 

«• 1/ 

k 1 

f 

* fc o 




While 

i < 

2. do : 







Let t 

= min (x^. 

t k ) 

in 


k 



a . 

If t 

* t' then: 







(1) 

P +- 
n 

<t>(t. 

f) 

P n 

4> T (t, 

t ’ ) + D* ( t , 

f) 


(2) 


$(t. 

f) 

9s 






8y 



9Y 





(3) 

91 

0(t, 

f) 

91 

+ 0(t, 

f ) 4 > a (t, 

t’) 



3B 



8B 


1 J 



(4) 

P «- 
u 

4>(t, 

f) 

P 

u 

0 T (t, 

t') 



+ <&(t, t') 

+ 0( t , t ’ ) 

+ 0( t , t * ) 
+ D(t , t * ) 


P uB ® T<t - t ' ) 

P uB <=■> 

d B 0 T (t , t*) 

- D* ( t , f) 
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(5) 


P uB * fc, > P uB * 1 ^' *=•> 

+ ©( t , f) d 0 <i>g(t, t’) + D sf3 (t, f) 

(6) dp 4 - 4>p(t , f) dp4>*(t, t') + dp(t, f) 

(7) <t> 4 - 4>(t, t*) 4> 

(8) d 4 - 4>(t, t’) d+ T (t, t*) + d(t , t*) 

b. If t = then: 

Let 

K = P G? (g. p gJ + ct 2 ) 
nk\knk k/ 

M = I - KG k 
in 

(1) P •»- MP M T + K a 2 K T 

n n k 

(2) ^ 4 M ^ - KI\ 

9y 9 y K 


(3) 

9s 

4- M ^ 


9(3 

93 

(4) 

P 

u 

4- MP M T 
u 

(5) 

P u(3 

4- MP „ 

u3 

(6) 

k 4- 
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c . 


If t = T i 

then: 


Let ** = 
3Y 

al_ • j 

J * * * ! 

Is 1 

QJ 

-< 

1 

3i = 

9(3 

Si_ S S 

93 1 ' * 1 

35 

9(3 

n (3 


in 

(1) P n < x i> *■ P n 

<■> w * E (%) (#;)’ 

( 3 ) 3 b <V * E (%) (337) 

(4) P u (T.) - P u 

(5) P^) «- P n (^i) + P Y< T i> + P f3 (x i ) + P u (T i ) 

(6) <t>(T i , x i _ 1 ) «- 4> 

(7) d(T., T i _ 1 ) 4- d 

(8) Reset d to zero and 4* to the N x N 

identity matrix 

(9) i i + 1 
d. t* «• t 


2-41 


0450 




SECTION 3 — SPIN— STAB ILIZED SPACECRAFT — MODELS 


As Section 2 discussed estimation filtering in a general 
sense, this section provides specific attitude and sensor 
models for spinning spacecraft. These models provide the 
basis for the construction of the state transition, measure- 
ment partial derivative, and random excitation matrices used 

in Section 2. 

3 . 1 ATTITUDE GENERATION 

For spin-stablized spacecraft , the attitude is determined by 
the direction of the spin axis. The nominal spin axis is 
assumed to be fixed in inertial space and is specified by 
its right ascension a n and its declination § n as shown in 
Figure 3-1. The actual spin axis may drift from this iner- 
tially fixed nominal position. ADEAS provides a model for 
sinusoidal variations of the actual spin-axis about the nom- 
inal axis: 

a* ( t) = a n + a a sin K (t - t Q ) + t> a ] (3.1-1, 

6 * ( t ) = S n + a § sin [£ (t - t Q ) + b g ] (3.1-2 

where the amplitudes a a and a g are constant dynamic param- 
eters, while the frequency C and the phases b a and b & are 
assumed to be exactly known. The spin-axis attitude defined 
by a* and 5* is used in all sensor model computations. 
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Figure 3-1. Right Ascension and Declination of the Nominal 
Spin Axis 

3 . 2 DYNAMIC ERROR MODEL 

If we define a state vector x by 


x 



(3.2-1) 


then differentiating Equations (3.1-1) and (3.1-2) gives the 
dynamic model for spin-stablized spacecraft in the form of 
Equation (2.1-1): 



'a a C cos [ C ( t - t Q ) + b a ]' 


1 

4J 

1 

x(t) = 

C cos [C(t - t Q ) + bg] 
0 

+ 

u 5 (t) 

0 


- 0 - 


- G - 
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where u and u, are dynamic noise processes that induce 

CL o 

additional attitude variations. 


From Equation (3.2-2), the state transition matrix, as de 
fined by Equation (2.1—5) , is then 


4>(t, t’) 


l o <J> o 

CL 

oio <j > 5 

0 0 1 0 

0 0 0 1 


(3 . 2-3a) 


with 


<{> a = sin [ C C t - t Q ) + b a ] - sin CC (t ’ - t Q ) + b a ] (3.2 3b) 


and 


4>- = sin [ C ( t - t n ) + b s ] - sin [C(f - t Q ) + b g ] (3.2 3c) 


The inverse of the state transition matrix is 


4> _1 (t, t') = 


1 

0 

0 

0 


0 

1 

0 

0 


-<t>. 


1 

0 


0 

-♦< 

0 

1 


(3.2-4) 


The a priori covariance of the attitude error parameters is 

p = diag ( ^ (3.2-5) 

^Aattitude \ a 5 / 
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Defining »J» as in Equation (2.1-6) the random excitation 
matrix is 


-Q a 0 0 

0 Q s 0 

0 0 0 

.0 0 0 

where 

E [“a (t> ^ (fc,) ] = Q a 6(t " t ’ ) 

E u $ (t) u£(f)J = Q 5 $(t - f) 

***** % 

The partitioning of d in Equation (3.2-6) is not the same as 
the partitioning in Equation (2.2-13); the two partitionings 
are related by row and column interchanges depending on the 
selection of dynamic solve-for and consider parameters. 

3.3 SENSOR MODELS 

The spin-axis sensors modeled by ADEAS are 

• IR horizon sensor 

• V-slit star scanner 

• V-slit Sun sensor 

For a spin-stabilized spacecraft, the primary attitude in- 
formation is given by the spin-axis direction U, which is 


01 

0 

0 

0J 


At 


(3.2-6) 


(3.2-7) 


d(t, t') = E 


vj/(t, t') y (t, t 


■>]■ 
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expressible in terms of its right ascension a and declina 
tion 6 in GCI coordinates as 


U 


I 


cos 5 cos a 

cos 6 sin a 

sin $ 


(3.3-1) 


The sensed data can be expressed in terms of a reference 
vector R, which is the unit vector in the direction of a 
sensed object. 

Since the spin-angle itself is not of interest /^the direc 
tion of the projection of the reference vector R in- the 
spin-plane conveys no attitude information. Meaningful 
measurements of any reference vector R give the following 
form of the measurement function h(x, p*) of Equa- 
tion (2.1-13): 


h(x, p 1 ) = h(C r , p*) (3.3-2) 

where h is a sensor-specific function 

and C is the projection of the unit reference vector R onto 
the spacecraft spin axis, or cosine of the angle between R 

AT 

and U ; 

C = R 1 cos 6 cos a + R* cos 6 sin a + R sin 5 (3.3-3) 

r x Y 

where the superscript I denotes GCI coordinates. 

For all spin-stablized spacecraft sensors, the functional 
dependence of the measurement on the spacecraft attitude, 
parameterized by the right ascension and declination of the 
spin axis, is through the projection C r - Since the attitude 
errors are always a subset of the solved-for vector, the 
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partial derivatives of the measurements with respect to the 
attitude are computed as 


3v* 2h_ 8C r 

3a 3C r 6a 


lh 

3C 




cos 5 sin a + R„^ cos 5 cos a (3.3-4a) 


3y*_ lh_ 

35 " 3C r 35 


( 3 . 3-4b) 


lh 

3C 


- -I 
r - 


sin 5 cos a - R sin 5 sin a + R cos 5 


The formulation of h and its partial derivatives with re- 
spect to the projection C r is thus a key to the sensor- 
related error analysis computations. The expressions for 
h(C r , p') are given below for each sensor as well as a com- 
plete list of all partial derivatives for each sensor. 

3.3.1 IR HORIZON SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The measurement model for the IR horizon sensor is shown in 
Figure 3-2. This model describes both the conical scan IR 
sensor, which has a fixed scan cone angle y, and the pano- 
ramic attitude sensor (PAS), which has a varying scan cone 
angle. Since the scan cone angle of a PAS is constant for 
each complete revolution of the spacecraft, the analysis can 
be carried out for a fixed scan cone angle. The analyst 
defines the parameters describing how the scan cone angle 
changes per revolution. If a zero incremental change is 
specified, the PAS will be identical to the conical scan IR 
horizon sensor. If the increment is not zero, the scan cone 
angle will increase each revolution until it reaches a user- 
defined upper limit and then decrease by the same incremental 
amount until it reaches the user-defined lower limit. This 
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Figure 3-2. IR Horizon Sensor Geometry/ Sun Indexed 
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> 10 / ‘ 8 ] 


is the only ADEAS sensor model that uses the projections 
onto the spacecraft spin axis of two external reference vec- 
tors, i.e., the Sun vector and the Earth vector. The meas- 
urements of the horizon sensor are the Sun to Earth-in 
azimuth. Sun to Earth-out azimuth. Earth width, and Earth 
azimuth with respect to the Sun. The measurements are cal- 
culated as 


y* = 


Earth width = h(Cg, C £ ) 


= 0 = 2 COS -1 j g - Q S - . P.- CQg.n CQS y | 


) 


sm n sm y 


f 


(3.3-5) 


y* = Earth azimuth = h(C_, CL.) 

O £* 


= 4> E = cos 


-1 I cos vjf - cos n cos ft 
) sin n sin 13 


(3.3-6) 


y* = Earth-in = h(C R , C F ) = A IN = 4^ - 1 n (3.3-7) 


Y* = Earth-out = h(Cg, C E ) = A QUT = ^ \ 0 (3.3-8) 


where Cg = projection of the unit Sun vector S^onto the 
spacecraft spin axis = cos (3 = §1 • 01' 

Cg = projection of the unit Earth vector g 1 onto the 
spacecraft spin axis = cos n = • U 1 

13 = Sun angle, angle between the spin axis and 
the Sun vector 

n = Earth angle, angle between the spin axis and 
the Earth vector 

cos = projection of the Sun vector on the Earth 
vector 

- 


p 

r 


= Earth angular 
-1 Ttl + 


radius 


= sm 
= ' Earth 


[ ( r + h t )/R] 
radius (km) 


h{- = IR tangent height (km) 
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r = magnitude of the spacecraft position vector 
(km) 

Y = sensor scan cone angle (deg) 

Equation (3.3-6) does not uniquely determine $ E ; the sign 
of the angle is not determined. To resolve this ambiguity 
the sign of a vector product is used. Assume that the cos 
function returns an angle in the range 0 < <I> E < 180. 

If (U 1 x S 1 ) • E 1 > 0 then Equation (3.3-6) gives the cor- 

rect range. 


AT 

If (U X 


) • E^ < 0 then use the following equation. 


<&„ = 3 60 - cos 

£j 


-1 j cos u> - cos n cos [3 I 
| sin n sin (3 ) 


(3.3-9) 


The measurement parameters, which the user may select as 
either solved for or considered, for the IR horizon sensor 
are 

1. The scan cone angle (deg) 

2. The Earth angular radius (deg) 

3. The distance from the Earth to the spacecraft (km) 

4. The IR tangent height (km) 

5. The Earth-in bias (deg) 

6. The Earth-in scale factor 

7. The amplitude of the Earth-in periodic measurement 
error (deg) 

8. The Earth-out bias (deg) 

9. The Earth-out scale factor 

10. The amplitude of the Earth-out periodic measurement 
error (deg) 
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11. The Earth width bias (deg) 

12. The Earth width scale factor 

13. The amplitude of the Earth width periodic measure- 
ment error (deg) 

14. The Earth azimuth bias (deg) 

15. The Earth azimuth scale factor 

16. The amplitude of the Earth azimuth periodic meas- 
urement error (deg) 

In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The initial sensor scan cone angle (deg) 

2. The incremental scan cone angle (deg) 

3. The minimum scan cone angle (deg) 

4. The maximum scan cone angle (deg) 

5. The IR tangent height (km) 

6. The frequency and phase angle of the Earth-in peri- 
odic measurement error 

7. The frequency and phase angle of the Earth-out 
periodic measurement error 

8. The frequency and phase angle of the Earth width 
periodic measurement error 

9. The frequency and phase angle of the Earth azimuth 
periodic measurement error 

When scheduling the IR horizon sensor the user may use one 
or two of the sensor outputs for the error analysis computa- 
tions. Since all outputs have units of degrees, the user- 
supplied value of the sensor white noise standard deviation 
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for each measurement has units of degrees. The errors in the 
Earth-in and Earth-out measurements are assumed to have equal 
standard deviations c A and to be uncorrelated with one 
another. Then, it follows from Equations (3.3-7) and (3.3-8) 
that the errors in the Earth— width and Earth azimuth measure- 
ments are uncorrelated, and that their standard deviations 
are 





'n 


= sTl < 7 , 


(3.3-10) 


The correlations between the pairs (ft, A^) • • ^OUT^ ' 

(<fr E , A in ), and (<I> E , A qut ) are ignored if any such pair is 

selected. 

Partial of the Earth azimuth wrt the projection of the Sun 

/N 

vector S: 


9<I> cos n - sin n cot (3 cos 4>. 
E _ «, 

ac g - sin n sin (3 sin 4> E 


(3.3-11) 


Partial of the Earth azimuth wrt the projection of the Earth 
vector E: 


9<t> cos 13 - cot n sin (3 cos ^ 
E _ & 

3C e “ sin n sin (3 sin <t> E 


(3.3-12) 


Partial of the Earth azimuth wrt the sensor scan cone angle: 


97 


0 


(3.3-13) 
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Partial of the Earth azimuth wrt the Earth angular radius: 



(3.3-14) 


Partial of the Earth azimuth wrt the IR tangent height: 



(3.3-15) 


Partial of the Earth azimuth wrt the distance from the 
spacecraft to Earth: 


3 ^ 

3R 


0 


(3.3-16) 


Partial of the Earth azimuth wrt the bias b: 



(3.3-17) 


Partial of the Earth azimuth wrt the scale factor k: 


34 > 
E 

3k 


~ h(Cg , C E ) = cos 


-1 I cos ^ - cos n cos (3 
) sin n sin 3 


(3.3-18) 


Partial of the Earth azimuth wrt the amplitude of the peri- 
odic error: 


30 e 

= sm (cot + 14 /) (3.3-19) 
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Partial of the Earth width wrt the projection of the Sun 
vector S: 


an 

ac. 


= o 


( 3 . 3 - 20 ) 


Partial of the Earth width wrt the projection of the Earth 


vector E: 


an 

ac T 


= 2 


cos y - sin y cos 9 cot n 


sin y sin n sin 


n 


( 3 . 3 - 21 ) 


Partial of the Earth width wrt the sensor, scan cone angle: 


M _ o 

3Y 


cos y sin n cos 9 - sin*-Y cos n 


. n 

sin y sin H sin 2 


( 3 . 3 - 22 ) 


Partial of the Earth width wrt the Earth angular radius: 


an _ , 

3p 


sin p 


n 

sin y sin n sm j 


( 3 . 3 - 23 ) 


Partial of the Earth width wrt the IR tangent height: 


an 

3h 1 


an 

3 P 


r — i — I 

|_R COS pj 


( 3 . 3 - 24 ) 


Partial of the Earth width wrt the distance from the space- 
craft to Earth: 


an _ _ an tan o 
3R ap R 


( 3 . 3 - 25 ) 
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Partial of the Earth width wrt the bias b: 


an 

3b 


= 1 


(3.3-26) 


Partial of the Earth width wrt the scale factor k: 

= h(C_, C_> = 2 cos' 1 j CTS ^ * ! (3.3-27) 

3k S E | sin n sin Y ) 

Partial of the Earth width wrt the amplitude of the periodic 
error : 


an 

3a 


= sin (cot + ifr ) 


(3.3-28) 


Partial of the Earth-in angle wrt the projection of the Sun 

/N * 

vector S: 


3A 

ac 


IN 

s 


3 ^ 

ac s 


(3.3-29) 


Partial of the Earth-in angle wrt the projection of the 

/N 

Earth vector E: 


3A 

ac 


IN 

E 


8c e 


1 an_ 

2 aC E 


(3.3-30) 


Partial of the Earth-in angle wrt the sensor scan cone angle: 


9a in _ _ i an 
3y 2 3y 


(3.3-31) 
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Partial of the Earth-in angle wrt the Earth angular radius: 


9a in 1 M 

3 p 2 3 p 


(3.3-32) 


Partial of the Earth-in angle wrt the IR tangent height: 


8A IN 1 3Q_ 

3h fc " 2 3h t 


(3.3-33) 


Partial of the Earth-in angle wrt the distance from the 
spacecraft to Earth: 


9A IN _ 1 3fl 

3R " 2 3R 


Partial of the Earth-in angle wrt the bias b: 


8a in 

3b 


1 


(3.3-34) 


(3.3-35) 


Partial of the Earth-in angle wrt the scale factor k: 


3A 


3k 


IN _ 


h < c s' c E ) 



(3.3-36) 


Partial of the Earth-in angle wrt the amplitude of periodic 
error : 



3a 


sin (ut + <J/) 


(3.3-37) 
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Partial of the Earth-out angle wrt the projection of the Sun 

/N 

vector S: 


9a out 

8c s 


3^ 

9C S 


(3.3-38) 


Partial of the Earth-out angle wrt the projection of the 
Earth vector E: 


9a out ^e i an 
ac E = ac E + 2 ac E 


(3.3-39) 


Partial of the Earth-out angle wrt the sensor scan cone 
angle : 


3A 

out i an 

3Y 2 3 y 


(3.3-40) 


Partial of the Earth-out angle wrt the Earth angular radius: 


9a out _ i an 

3p 2 3 p 


(3.3-41) 


Partial of the Earth-out angle wrt the IR tangent height: 


9a 

OUT _ 1 an_ 
3h. 2 3h fc 


(3.3-42) 


Partial of the Earth-out angle wrt the distance from the 
spacecraft to Earth: 


9a out _ i an 

3R 2 3R 


(3.3-43) 
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Partial of the Earth-out angle wrt the bias b: 


3A OUT 

9b 


(3.3-44) 


Partial of the Earth-out angle wrt the scale factor k: 


Partial of 
error : 


9A 


OUT 


9k 


« h(Cg, C E ) = 4> p + ^ £2 


E 2 


the Earth-out angle wrt the amplitude of 


(3.3-45) 

periodic 


- 7 ™ = sin (cot + Y) (3.3-46) 

d a 

3.3.2 V-SLIT STAR SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The V-slit star sensor contains two slits, one parallel to 
the spin axis, the other at an oblique angle, as in Fig- 
ure 3-3. The raw measurement is the time At between the 
crossing of each slit by the star image focused by an op- 
tical system. ADEAS uses as its preprocessed measurement y* 
the rotation angle of the spacecraft between the star cross- 
ings, given as coAt where co is the angular rotation rate of 
the spacecraft. This measurement is a function of the cosine 

A J 

of the star angle a between the star vector S and the 
spin vector as 

y* E M = rotation angle = h(Cg) (3.3-47) 

= sin -1 [tan (y - °) tan f] + A 
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Figure 3-3. V-Slit Star Scanner Spherical Geometry 


where Cg = projection of the unit Sun vector S^onto the 
° spacecraft spin axis = cos a = S 1 • U 

a = star angle 

Y = mounting angle of the optical axis of the star 
scanner 

E = tilt of oblique slit 

A = separation of center of slits, the angle between 
the two lines perpendicular to the spin axis from 
that axis to the intersections of the slits with 
the plane scanned by the optical axis 

Though this expression is exact only when the mounting angle 
is 90 degrees and there is no rotation of the star scanner 
about its optical axis, it is a good approximation when the 
lengths of the slits are small, typically 10 degrees. 

To determine the sensitivity of the measurement with respect 
to a rotation of the slits about the optical axis, a differ- 
ent measurement model is used. 

y* = rotation angle = h(C g ) = sin -1 [tan (y - <*) t an (2 + 0)] 

- sin -1 [tan (y - °) tan 0] (3.3-48) 

B A - B 

+ cos 0 + (1 - tan © tan E) cos 0 

where 0 = rotation of the slits about the optical axis 

B = angular distance between the optical axis and 
the vertical slit 

The last two terms of this model are based on plane geometry 
(Figure 3-4 ) . 

The geometrical limitation on the visibility of the V-slit 
star sensor requires that the star angle be between the min- 
imum and maximum values of 


a . 
min 



(3.3-49) 
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(3.3-50) 


1 

a = y + - A. 
max 2 

where X is the length of the vertical slit. 

The measurement parameters, which the user may select as 
either solved for or considered, for the V-slit star sensor 

are 

1. The measurement bias (deg) 

2. The measurement scale factor 

3. The amplitude of the measurement periodic error 
(deg) 

4. The separation of the center of the slits (deg) 

5. The tilt of the oblique slit (deg) 

6. The rotation of the slits about the optical axis 
(deg) 

Even though the rotation of the slits about the optical axis 
is always nominally 0 degrees, it may have a specified un- 
certainty and hence be designated as an error parameter. 

In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The frequency and phase angle of the periodic meas- 
urement error 

2. The separation of the center of the slits (deg) 

3. The tilt of the oblique slit (deg) 

4. The angular distance between the optical axis and 
the vertical slit (deg) 

5. The length of the vertical slit (deg) 
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When scheduling the V-slit star scanner, the analyst does 
not have a choice of measurements to use, since there is 
only one. Since the rotation angle has units of degrees, 
the sensor white noise standard deviation has units of 
degrees . 

Partial of the rotation angle M wrt the projection of the 

/\ 

star vector S: 


3M = 1 3M 

3C sin a 3y 
s 


(3.3-51) 


Partial of the rotation angle M wrt the separation of the 
center of the slits: 


3M 

3A 


(3.3-52) 


Partial of the rotation angle M wrt the tilt of the slit: 

M m tan (y - a) 3 

3E cos 2 E cos (M - A) 

Partial of the rotation angle M wrt the mounting angle of 
the optical axis: 


3M 

3y 


tan E 

2 

cos (y - o ) cos (M-A) 


(3.3-54) 


Partial of the rotation angle M wrt rotation of the slits 
about the optical axis, from Equation (3.3-43): 


QQ = (A - B) tan E + tan (y - a) tan 2 E (3.3-55) 


3-22 


0450 



in the limit 0 -*■ 0 and y - o « 1. 

Partial of the rotation angle M wrt the bias b: 

9M _ , (3.3-56) 

3b " 

Partial of the rotation angle M wrt the scale factor k: 

a h(C g ) = A + si n_1 [tan (y - <*) tan E] (3.3-57) 

Partial of the rotation angle M wrt the amplitude of the 
periodic error: 


~ = sin (cot + Y) (3.3 58) 

9a 

3.3.3 V-SLIT SUN SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The V-slit Sun sensor contains two slits of equal length, 
one parallel to the spin axis, the other at an oblique 
angle, as in Figure 3-5, from which is determined the time 
At between Sun crossings. ADEAS uses as its preprocessed 
measurement y* the rotation angle of the spacecraft between 
Sun crossings given as coAt where to is the angular rotation 
rate of the spacecraft. This measurement is a function of 
the cosine of the Sun angle (3 between the Sun vector S and 
the spin vector. 

The rotation angle is given by 

M = sin -1 (cot 3 tan E) - sin -1 (cot 3 tan Y) + A (3.3-59) 
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in the general case illustrated in Figure 3-5. ADEAS only 
models the special case of t s 0 , in which case the measure- 
ment equation is 

y * = M = rotation angle = h(C g ) = sin -1 (cot 13 tan Z) + A 

(3.3-60) 

where C s = projection of the unit Sun vector S 1 onto the 
spacecraft spin axis = cos 13 

Z = tilt of the oblique slit 

A = angle between the lines formed by the intersec- 
tion of each slit with the x-y plane, the coor- 
dinate frame being the spacecraft's principal 
axes, with the spin vector parallel to the z-axis 

This is the same as the V-slit star sensor model. Equa- 
tion (3.3-47), with y = ir/2 and a = (3. 

If the two slits have the same lengths, -a-s is assumed in 
ADEAS, the slit with the greatest tilt (I) will define the 
maximum and minimum (3 that .can be measured. Thus, the geo- 
metrical limitation on the visibility of the V-slit Sun sen- 
sor requires that the Sun angle (3 be between the minimum and 

maximum values of 

a =90 deg - sin -1 (sin (X/2) cos Z) deg (3.3-61) 

min 

o =90 deg + sin -1 (sin (X/2) cos E) deg (3.3-62) 

max 

where \ is the length of the slits. 

The measurement parameters, which the user may select as 
either solved for or considered, for the V-slit Sun sensor 

are 

1. The measurement bias (deg) 

2. The measurement scale factor 
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The amplitude of the measurement periodic error 
(deg) 

4. The angle A between the lines formed by the inter- 
section of each slit with the x-y plane (deg) 

5. The tilt of the oblique slit (deg) 

6. The tilt of the vertical slit (deg) 

Note that even though the tilt of the vertical slit is 
always nominally 0 degrees, it may have a specified uncer- 
tainty and hence be designated as an error parameter. 

In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The frequency and phase angle of the periodic 
measurement error 

2. The angle A between the lines formed by the inter- 
section of each slit with the x-y plane (deg) 

3. The tilt of the oblique slit (deg) 

4. The length of the slits (deg) 

When scheduling the V-slit Sun sensor, the analyst does not 
have a choice of measurements to use, since there is only 
one. Since the rotation angle has units of degrees, the 
sensor white noise standard deviation has units of degrees. 

Partial of the rotation angle M wrt the projection of the 
unit Sun vector S: 


9M _ tan Z 

8C s sin 3 (3 cos (M - A) 


(3.3-63) 
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Partial of the rotation angle M wrt the angle between lines 
of slit intersection with the xy plane: 

= 1 (3.3-64) 

9 A 

Partial of the rotation angle M wrt the tilt of the tilted 
slit : 


225 = cot_JB (3.3-65) 

8Z cos 2 E cos (M - A) 

Partial of the rotation angle M wrt the tilt of the vertical 
slit : 


(3.3-66) 
(3.3-67) 

Partial of the rotation angle M wrt the scale factor k: 

If ~ h(C ) = A + sin -1 (cot 13 tan E) (3.3-68) 

3k s 

Partial of the rotation angle M wrt the amplitude of the 
periodic error: 


= -cot (3 


Partial of the rotation angle M wrt the bias b: 


3M 

9b “ 


9M 

9a 


sin (cot + 


(3.3-69) 
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SECTION 4 - THREE-AXIS STABILIZED SPACECRAFT M ODELS 


As Section 2 discussed estimation filtering in a general 
sense, this section provides spe'cific attitude and sensor 
models for three-axis stabilized spacecraft. These models 
provide the basis for the construction of the state transi- 
tion, measurement partial derivative, and random excitation 
matrices used in Section 2. 

4.1 ATTITUDE GENERATION 

For three-axis stabilized spacecraft, the attitude is defined 
by the 3x3 orthogonal transformation matrix A B/R from some 
reference coordinate system (see Section 5.1) to the space- 
craft body coordinate system. The nominal spacecraft at- 
titude A’ , evolves over time based on the body components 

of the nominal spacecraft angular velocity u^ /R relative to 

the reference coordinate system (see Section 16.1 of Refer- 
ence 3) : 


A B/R^ " “B/R^ A B/R^ 


(4.1-1) 


where coJ, /Ti is the 3x3 antisymmetric matrix defined from 
B/R 

“b/r by: 


“b/r 


0 

-a) B/Rz 

“B/Ry 

W B/Rz 

0 

- W B/Rx 

-t0 B/Ry 

U B/Rx 

0 


(4.1-2) 


In addition, ADEAS models sinusoidal stabilization errors 
about the nominal attitude: 


A B/R^ “ ^e^ A B/R^ 


(4.1-3) 
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where 


♦e (t > 


I + 


reftn 2 

|e(t)| 2 


[1 - COS 


I e(t) | ] 


sin |e(t)| 

|e<t) | 

(4.1-4) 


and e is an antisymmetric matrix defined from a vector e 
similar to Equation (4.1-2). The variation vector e is 
given by 


e i = a i s * n ^i ^ " t Q ) + ^i^ i = 1, 2, 3 (4.1-5) 

where the amplitudes a^ the frequencies C.^, and the phases 

v. are all constant. The attitude matrix 
1 1 

A B/I = A B/R ^/I (4.1-6) 


is used for all sensor model computations, where the 
inertial-to-ref erence matrices A R/I for different reference 
frames are found in Section 5.1. 


For ADEAS, the nominal attitude A R ^, R (t) is specified by an 
initial attitude and a set of attitude rates. The initial 
attitude A ’ /tl (t n ) may be specified by Euler angles or a 
quaternion (see Reference 3, Section 12.1 and Appendix E) . 


The angular velocity Wg /R is specified between discrete times 
t^ in one of two ways. The first is as a set of constant 

The Euler angles for 


Euler angle rates, 4>, 9, y. 
t^_^ < t < t^ are given by 


4>(t) = 4>(t i _ 1 ) + ^(t - t i _ 1 ) 

0(t) = 0(t i _ 1 > + e.<t - t._ 1 > 

v( t) = »♦»( t i _ 1 ) + <i» i (t - t i _ 1 ) 


(4.1-7) 
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If the initial attitude is specified by Euler angles and co£ /R 
by Euler angle rates, the same Euler axis sequence is used 
for both. If the initial attitude is specified by a 
quaterion and c3£ /R by Euler angle rates, the initial Euler 
angles 4><t Q ), 0(t Q ), and y(t Q ) are computed from the initial 
attitude matrix A£ /R (t Q ) using the equations in Table 4-1. 

The correct quadrants for 4> and are determined by the fact 
that the signs of the numerator and denominator of the argu- 
ment of tan" 1 are the signs of the sine and cosine, respec- 
tively, of the angle. If both the numerator and denominator 
of the expression for tan 4> are zero, <j> is set to zero. 

The second means of specifying w£ /R between the discrete times 

t. is as a constant vector of components in the body frame. 

When oi' _is constant, the solution to Equation (4.1-1) is 
B/ R 

*4, R (t > * *A (t - t i-l ) A B/R <4 - 1 ‘ 8> 


where wA /t> is constant for ti_i < t < ti and 
B/ K 


+A (t ' t i-l > 


1 2 


= I + _ 


IfiZEi 




[1 - COS I^b/rI (t - 


fc i-l> ] 


(4.1-9) 


^B/R 

I“b/r! 


sin I w b/R I (t " t i-l ) 


4.2 DYNAMIC ERROR MODEL 

Let A 0/r be the true spacecraft attitude and A* /R be an es- 
timate of the attitude. For all internal processing, ADEAS 
represents the error in A* /R by the first order error vec- 
tor A0 . The components of this vector represent the small 
rotations needed about each of the spacecraft body axes to 
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of Euler Angles From Attitude 
3 ) 


4-4 



Table 4-1. Computation of Euler Angles From Attitude 
Matrix (2 of 3) 


i-i-k 

3-1-2 


3 - 2-1 


1 - 2-1 


1-3-1 


EULER ANGLES 


<j> = tan 1 

[ A 21 /A 22] 

0 = sin 1 (A 23 > 

-1 

"A^ 7 sin 4> + A 31 cos <t> 

y = tan 

A 12 sin <t> + A i;l cos <t> 


4 > 

0 


tan 1 jA 12 /A ii] 


sin 1 ^“ A i 3 ^ 


y = tan 


-1 


A 31 sin 4> - a 32 cos 4> 


L - A 2 i sin 4> + A 22 cos <t> 


4 > = tan -1 

[ A 12 / A 13^ 

0 = cos 1 ( a i;l ) 

-1 

sin <t> - A ^ 2 cos 4> 

v|f = tan 

A 23 sin 4 > + A 22 cos 4> 

4 > = tan - ”* - 

[ A 13 /A 12] 

0 = cos 1 < 

-1 

P-A 22 sin $ + A 23 COS ^ 

t)i = tan 

[-A 32 Sin a + A 33 cos 4* 
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Table 4-1 


Matrix (3 


2-1-2 <t> = 

9 = 

Y = 

2- 3-2 4> = 

e = 

Y = 

3- 1-3 4> = 

e = 

Y = 

3-2-3 4> = 

© = 

Y = 
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align them with the estimated body axes (see Figure 4-1) . 

The true and estimated attitudes are then related by 

a |/e « < T + A§ > a b/r (4 - 2 - 1) 

where the antisymmetric matrix A0 is defined similarly to 
Equation (4.1-2). 

4.2.1 ATTITUDE ERROR PARAMETERIZATION 

While the error vector A© is useful for internal computa- 
tions/ ADEAS provides two additional attitude error param- 
eterizations for convenience of input and output: Euler 

angle errors and quaternion errors. 

Euler angle errors represent the difference in the Euler 
angles between the Euler sequence rotating the reference 
coordinate system to the true spacecraft body system and the 
Euler sequence rotating the reference system to the estimated 
body system. The explicit forms for the attitude matrix in 
terms of the 12 possible Euler sequences are given in 
Table E-l of Reference 3. The relationship between the 
attitude error vector and the Euler angle errors is analogous 
to the relationship between the angular velocity vector and 
the Euler angle rates. This analogy holds because the atti- 
tude errors and the actual angular motion in an infinitesimal 
time At are both infinitesimal. Therefore, the transforma- 
tions between the components of the attitude error vector A0 
and the Euler angle errors A<t>, A0, Ay can be obtained from 
the well known relations relating Euler angle rates and the 
angular velocity components. The following table of trans- 
formations (Table 4-2) for 12 different Euler sequences is 
adapted from Reference 3. 


4-7 


0450 













The three columns of the inverse transformation matrix rep- 
resent the axes of Euler rotations expressed in the orthog- 
onal body axes. It should be noted that the adjacent columns 
are orthogonal. Similarly, the adjacent rows of the trans- 
formation matrix are orthogonal. The Euler angles <t>, 0, and 
\|r are derived from the attitude matrix by the transfor- 

mations in Table 4-1. 


The quaternion errors represent the difference in the qua- 
ternion describing the rotation from the reference coordi- 
nate system to the true spacecraft body system and the 
quaternion describing the rotation from the reference system 
to the estimated body system (Section 12.1 of Reference 3). 
As in the case of Euler angle errors, the relationship to 
the attitude error vector is analogous to the relationship 
between the angular velocity and the quaternion rates. 
Multiplying the dynamic equations of motion for quaternions 
by At gives (see Section 16.1 of Reference 3) 


-Aq’ 


" q 4 

- q 3 

q 2~ 

Aq 2 

1 

q 3 

q 4 

" q l 

Aq 3 

“ 2 

- q 2 

q l 

q 4 

- Aq 4_ 


r q i 

- q 2 

“ q 3 



(4.2-2) 


where the quaternion q = 
the attitude matrix A B/R 


(q 1# q 2' q 3' q 4^ is derived from 
computed as in Section 4.1. 


As stated previously, the error vector A0 is used for inter- 
nal computations. Thus, the a priori covariance of the at- 
titude error parameters is always given by 



P Aattitude - E < 4 - 2 ' 3a > 
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The assumption is made in ADEAS, however, that the initial 
uncertainties in the parameters used for input and output 
are uncorrelated. Thus, the form of the a priori covariance 
of the attitude error parameters depends on which parameteri- 
zation is used for input and output. 

If the attitude error vector is used for input and output, 
then 


P Aattitude - diag (V V 4 3 j 


( 4 . 2-3b) 


If Euler angle errors are used for input and output, then 


P Aattitude B 


diag 


(<V a Q’ a £)] 


B 


( 4 . 2-3c) 


where B is the appropriate matrix from the last column of 
Table 4-2. 

The situation is more complicated if quaternion errors are 
used for input and output since the quaternion errors are 
not independent but must obey the normalization condition 


Aq = 0 ( 4 . 2-3d ) 

Thus, we assume that the quaternion errors are given by 


Aq = (I - q q^Sq (4 . 2-3e) 

where the components of Sq are assumed to be independent so 
that 


EtSqSq 1 ’] 


diag 
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With the Aq given by Equation (4.2-3e), Equation (4.2-3d) 
is satisfied automatically. It would be inconsistent with 
Equation (4.2-3d) to assume that E[AqAq ] had the form of 
the right side of Equation (4.2-3f). Now we find from Equa- 
tions (4.2-2), (4 . 2-3a) , (4.2-3e), and (4.2-3f) that 


P 


Aattitude 



diag 






B 


q 


(4 . 2-3g) 


if quaternion errors are used for input and output, where 
B is the 4x3 matrix appearing in Equation (4.2-2). 

q 

4.2.2 ATTITUDE ERROR PROPAGATION 


The true attitude A fi/I relative to inertial space, given by 
Equation (4.1-6), evolves according to 


A B/I^ ” ““b/I^ A B/I^ (4.2 4) 


where 


(d 


B/I 


" “b/R + “r/I 


(4.2-5) 


is the angular velocity of the spacecraft relative to iner- 
tial space, with w R/I being the angular velocity of the ref- 
erence coordinate system.- All angular velocity components 
in Equation (4.2-5) are in body coordinates. Similarly, the 
estimated attitude 


A B/I " A B/R ^/I 


(4.2-6) 
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evolves according to 


A*/i(t) = -S* /T (t) A S/I (t) 


B/I 


B/I 


(4.2-7) 


where w* is the column vector of angular velocity compo- 
B/I 

nents in the body coordinate system derived using the gyro- 
scope measurement model (see Section 4.3.1). 

From Equation (4.2-1) and using Equations (4.1-6), (4.2 1) , 

and (4.2-4) through (4.2-7) we have 


dt 


dt ( a b/r a b/r) dt ( a b/i a b/i) 


= a b/i a b/i + a b/i a b/i 


(4.2-8) 


= _ “b/i a b/i a b/i + a b/i a b/i u b/i 


= - “b/I A0 + A0OJ B/I “b/I + W B/I 


Let 


A(0 


B/I 


= ( 0 * 


B/I 


- CO 


B/I 


(4.2-9) 


and assume that AtOg^ is small. Then, to first order 


dt 


A© = 


-td B/I A0 + A0 “B/I A “b/I 


(4.2-10a) 


Or, in vector form 


A0( t) = -“ B /I (t) A0(t) “ A “B/I (t) (4.2-10b) 
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According to the gyroscope model of Section 4.3.1 


Ao) B/I (t) = Ab( t ) + Q(t) Ak - “ B/I <t) Ae - u Q (t) (4.2-11) 

where Ab is the gyro bias error (a first-order Markov proc- 
ess), Ak is the gyro scale factor error, Ae is the gyro 
alignment error vector, Uq is a white noise process and 

fl( t) E diag [w 0/I (t)] (4.2-12) 

The attitude error A0 is nominally zero and thus can be in- 
cluded with Ab, Ak, and Ae in a composite state error vector 


A0 

Ab 



Ae 


(4.2-13) 


Combining Equations (4.2-10), (4.2-11), and (4.3-22) then 
gives the following state error equation: 




_0J B/I 

(t) -I 

-Q(t) 

S B/I < fc > _ 


1 

1 c I 
CD 
✓**>. 
rr 

v/ 

1 

• 

Ax(t) 



0 

— I/T 

0 

0 

Ax(t) + 

u b (t) 



0 

0 

0 

0 


0 



0 

0 

0 

0 


- 0 _ 







(4.2-14) 

where 

u fa (t) is 

a white noise 

process 

• 
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4.2.3 TRANSITION MATRIX COMPUTATION 

Integrating Equation (4.2-14) then gives for the state tran- 
sition matrix, as defined by Equation (2.1-5) 


<Mt, f) = 


n> ee ( t, t-> 

+9b< t - t-> 

1 

1 

1 


I - * ee <t. f) 

0 

*bb< fc ' *•> 

1 

1 

1 

1 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

1 

1 

1 

0 

I 


(4.2-15) 


where 


♦ee<t. c> - - S B/i (t) W’ t ' ) 


(4.2-16) 


4> 0b <t, f) = -S R/T (t) 4> 0b <t, f) - ie 


— ( t— t ' ) /t 


B/I 


(4.2-17) 


♦ 0k (t, f) « -S R/T (t) 4> 0k (t , f) - Q(t) (4.2-18) 


'B/I 


4>bb< fc ' f) - ie 


-(t-t* )/T 


(4.2-19) 


This partitioning of the transition matrix is different from 
the partitioning of Equation (2.2-ld); the two partitionings 
are related by row and column interchanges, depending on the 
selection of dynamic solve-for and consider parameters. 

Now Equation (4.2-16) has an identical form as Equa- 
tion (4.2-4) for the attitude A B/I . Thus, 4> 00 must also 
act as a transition matrix for • 


(t) = Wt, f) A fi/I (f) 


a b/i vw = ^ee 


(4.2-20) 
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or 


i 


^(t, f) = A B/I (t) A B/I ( t ' ) (4.2-21) 


Given a solution for 4 >qq, Equations (4.2-17) and (4.2-18) 
can be integrated to give 


^t ’ 


4> eb <t, fc, > = - / W*' fc "> e_<t " t >/T dt " < 4 - 2 " 22 > 


^(t, t 


■> - -/ 


♦qqU, t") £l(t") dt" (4.2-23) 


Substituting Equation (4.2-21) into these gives 


^©b^' = a b/I^ 


Ag /T ( t " ) e - ^ t " -t ' > /T df 


V 


(4.2-24) 




V 


= “ A B/I ( t) Aj(t. f) 


<j> ek (t, f) = -A B/I (t) / A£ /;[ (t") Q(t") df 

* / t ' 


- - A B/I (t) A k (t ' fc ’> 


(4.2-25) 
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where 


A fa (t, t 


y t ' 


A u/T (t”) e _<t ’ , " t ' )/T dt" (4.2-26) 
B/ 1 


A k (t, t 


•>■/ 

y t ’ 


Q(t M ) A B/I (t") dt" (4.2-27) 


Over a small interval At such that l“ B/I l <<: 1 ' the tran 
sition matrix can be approximated, to first order m 


“ B /il At ' by 


4> ee (t, t - At) ~ I - S B/I (t) At 


(4.2-28) 


Assuming also that At/x « 1, then, to first order. 


e -At/x ~ x _ At/X 


(4.2-29) 


Using Equations (4.2-28) and (4.2-29) in Equation (4.2-22) 
gives 

♦ 0b (t, t - At, * - f 1 1 - * **»( l - ‘ x*^) dt ' 

t-At 

= - f [(i - ** - x f ** r )- =B/I <t>(t - **>] dt ‘ (4.2-30) 


t-At l 

-IAt ♦ \ I ♦ S B/I (t) ] “ 2 


Note that this result is still first order in |« B/I I At and 
At/x even though it is second order in At. Using Equa- 
tions (4.2-28) and (4.2-29) again gives 


4> eb (t, t - At) 


- \ At[le- 4t/T ♦ * 0e (t, t - At)] (4.2-31) 


4-17 


0450 



Or, by Equation (4.2-21) 


W fc . t - 4t> * - \ At 


. —it/ T 


♦ A i/I< t - 4t) 


■] 


(4.2-32) 


Then, from Equation (4.2-24) 


A b (t, t - At) « \ At i 


A B/I^ fc ^ e + A B/I^ fc 


B/I 


(4.2-33) 


Substituting Equation (4.2-28) into Equation (4.2-23), and 
taking 0 to be approximately constant over the interval At, 
gives 

<t> ek (t, t - At) ~ 


where we have written Q a „ 0 (t) to emphasize that this is 
really the average Q over the time interval from t - At to t. 

Using Equations (4.2-28) and (4.2-21) then gives 


t 

t-At 




[I - S B/I (t)(t - t")] G ave <t) dt" 




c 


(4.2-34) 


lit - £ S B/I (t) it 2 | n a „„(t) 


1 


ave 


4. ek (t, t - it) 


- in - +ee (t - t - 4t)1 Q ave< k > 4t 

" 2 Afl/I* 4 * [ A B/I^ 

♦ Ag/I < t - ^)] a ave (t> it 


(4.2-35) 
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Then, from Equation (4.2 25) 


A k (t, t - At) ~ £ fl ave (t) At[A B/I (t) + A B/I (t “ At>1 (4 ' 2 36) 

In many cases, it is necessary to propagate the covariance 
over an interval t - t* that is large enough so that either 
or both of Equations (4.2-28) and (4.2-29) would be violated 
for At = t - f. In this case, the full interval is broken 

up into n equal steps of length: 

At = (t - t*)/n (4.2-37) 


where n = 1 + TRUNC{max[ | p ( t , t')|, (t - t')/r]/S} (4.2 38) 

5 is a user-specified tolerance, TRUNC is the truncation 
function, i.e., TRUNC(x) is the largest integer less than or 
equal to x, and p(t, f) is the net rotation of the space- 
craft with respect to inertial space between times f and t. 
If gyro biases are neither solved-for nor considered, 1/x 
is set to zero, so that it does not affect the number of 
steps. Equations (4.2-37) and (4.2-38) then guarantee that 

max( |w ave l' At, At/x) < 6 (4.2-39) 

where the average angular velocity to ave is related to 
the total rotation by 


p(t. f) = <0 ave (t - f) 


(4.2-40) 
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The vector p(t, t') is computed by using Equation (4.2-40) 
and the analogs of Equations (4.1-7) and (4.1-8) for rota- 
tions with respect to inertial space, which give 


^ 2 

A_ . _ ( t ) A*(f) - I + t - COS |p(t, t * ) | ] 

B/I tJ/1 „ /+. 4- • \ ^ 


|p(t, t’)| 

- P (t * sin |p(t, t')| 

I p(t, t ’ ) | 


(4.2-41) 


The trace of this equation gives 


p(t,t ’ ) | = cos 


-1 


1 

2 


tr 


rWi^ 


B/I 


(f 


>] -i 


(4.2-42) 


Note that |p(t,t')| represents the rotation through the 

smallest angle connecting the initial and final attitudes, 

rather than the actual rotation performed in a maneuver. 

For a complex maneuver, it may be that l w ave l At « 

I u I At > 6 . 

1 1 ave 

The computations used in ADEAS will be accurate if either a 
measurement or an output is processed between different arcs 
of a complex maneuver. 

The matrix fl (t) At, which is needed in Equation (4.2-36), 
is given, using Equation (4.2-12), by 


Q ave^ At = At] (4.2-43) 

where to_, T (t) means the average angular velocity from t - At 
B/ 1 

to t . 
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The nonzero elements of Q ave <t)At are numerically equal to 
three of the elements of the antisymmetric matrix « B/I <t)At, 

which can be computed from A B/I (t) and A B/I (t - At) by taking 
half of the antisymmetric part of Equation (4.2-41) with 
t' = t - At: 

2 [ A B/I (t " At)A B/I (t) ‘ A B/I (t) A B/I (t ' At> ] 

, 7T(t.t - &tl sin - At)| (4.2-44) 

|p(t,t - At)| 

«p<t.t - At) - « B/I (t)At 

The approximation is valid if |p(t,t - At) | is small 

enough so that |p(t,t - At) | « sin |p(t,t - At) | , 

which is guaranteed by the choice of 6 to be at least as good 

as other approximations made in the dynamic analysis. 

Now the integrals for A fa and A R in Equations (4.2-26) and 
(4.2-27) can be approximated by the sum of the corresponding 
values over each of the intervals At/ as given by Equa 
tions (4.2-33) and (4.2-36). Thus, letting = t* + iAt: 


n-1 

£ 

-iAt/x , /. e. \ 

e A b (t i+1' t i ) 


i = 0 


n-1 


-1 


e ~iAt/x 

A B/I (t i+1* 

-At/x / t . ) 

e + A B/I 


i = 0 


(4.2-45) 



n-1 
i = l 


A B/I ^ ^ + A B/I^ fc ^ e 


-nAt/x | 


A B/I (t i> e 


-iAt/xl 
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n-1 


A k (t, f) 


- Z V'i.l' •=!> 

i = 0 

(4.2-46) 

n-1 

= iZ °ave (‘i.P + A B/I (t i )I 


i = 0 


The transition matrices ^^(t, t’) and t’) 

can then be computed using Equations (4.2-24) and (4.2-25). 


Since the submatrix ^q q ( t , t‘) is seen from Equa- 
tion (4.2-21) to be orthogonal, the inverse of the full 
state transition matrix is given by 


♦ee (t ' fc,) 

*•> »•> «•> 

-♦Jett. fc ’> *’> 

i - ♦Jett. t*> 

(4.2-47) 

0 

»•> 

0 

o 

o 

! o 

i 

0 

0 

I 0 

0 

z 


where 

<t>b£(t, t* ) = ie (t-t,)/T (4.2-48) 
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4.2.4 RANDOM EXCITATION COVARIANCE MATRIX COMPUTATION 

Defining ^ as in Equation (2.1-6)/ then, based on Equa- 
tions (4.2-14) and (4.2-15), the random excitation matrix is 

d(t , t‘) = E[^(t, t’) ?’(t, t’>] 
t 

4>(t, t") Q u <t> T (t, t") dt M 

(4.2-49) 


d ee (t ' ^ 

1 

1 

1 

a © b (t. f) 

1 

» 

* 

o ! 

— f. 

0“ 

a eb (t ' t ’ > 
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1 

1 

1 

a bb <t, f) 
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1 

1 

I 
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o ! 

L- 

0 

0 
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0 

r“ 

i 

i 

o i 

0 


-L 



— !_ 

_ _ _ _ 

0 

t 
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0 

T" 

1 

1 

o ! 

0 _ 


where 

f> q 9 f) 

(4.2-50) 

•) o b <t>e b ( t - f>] dt- 
") Q b 4> bb (t, t") dt" (4.2-51) 

") Q b 4> bb (t, t") dt” (4.2-52) 


d ee (t, t‘) = / |Wt, 


©e' 


V 


+ W*' 


a eb (t 


,f>-/ 

•v 


+eb (t ' t ' 


a bb (t ' t ' ) - / W*' 
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with 


l 


E^u e (t) u^t')] = Q 0 5(t - f) (4.2-53) 

E [u b (t) ^<t')] = Q b 8 ( t - f) (4.2-54) 

Note that this partitioning of d is different from the par- 
titioning of Equation (2.2-13). The two partitionings are 
related by row and column interchanges, depending on the 
selection of dynamic solve-for and consider parameters. 

Substituting Equations (4.2-19), (4.2-21), and (4.2-24) into 
Equations (4.2-50), (4.2-51), and (4.2-52) gives 



(t, f) = A fi/I (t) d^t, t') Ag /Z ( t ) (4.2-55) 

(t, f) = -A B/I (t) a 0b (t , f) (4.2-56) 


J T[1 - e 2<t fc,)/T ] Q b , if 1/ x * 0 


d bb (t, f) = 


Q b (t - t') 


, if 1 /t = 0 


(4.2-57) 


where 


d£ e (t, t 



ij/i 


(4.2-58) 


+ A b (t, t") Q b A b (t, t")] dt" 


4-24 


0450 



Over a 
we can 
tions 

d ee (t/ 


l 0b (t ' 


t’) 


■/ 


A b (t ' 


) e -(t-t")/T Q dt « (4.2-59) 


small interval At, with |w B/I l At « 1 and At/x « 1 
use the first order approximations of Equa- 
4.2-28) and (4.2-30) in Equation (4.2-50) to obtain 


t - At) f |q@ + [Q@ “b/I^ “b/i^ Q e ] (t fc ^ 

*4-At ( 

♦ Q b (t - f) 2 + ^[Q b S B/I <« 

-s B /x< t > Q b I < t - t,)3 

-iv i - t >3 } at- 

= Q e At + ^[Q e w B/I (t) - w B/I (t) Q e ] At 2 

+ 3 Q b Afc3 + 8 [Q b “B/I (t) " S B/I (t) Q b ] At4 
- Q b 4X - 

Atj^Q e + 4> e0 (t, t - At) Q 0 4>| 0 (t, t - At) 

+ £ Q b (At) 2 + \ <fr eb <t, t ~ At > Q b *eb (t ' fc " At> ] 

(4.2-60) 
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Substituting Equations (4.2-21) and (4.2-24) into this gives 


d ee (t ' * - At> s 2 “ A B/i (t) [ A B/i (t > Q e A B /i (t> 

+ Ag /X (t - At) Q 0 A B/I (t - At) 

* 6 A B/I (t > «b A B/I (t)(At > 2 <4 ' 

. 1 A^(t, t - it) Q b A b (t, t - at)] Ag /r 

Then, from Equation (4.2-55) 

d ee (t ' t - at) * 1 at [A^ /r <t) q 0 A B/I (t) 

+ A B/i (t - 4t) Q a a b/i < - “> 

(4. 

+ i A Ll (t > Q b A B/I (t > <4t > 2 

+ \ A^(t, t - at) Q b A b (t, t - at)] 
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Using the approximations of Equations (4.2 29) and (4.2 30) 
in Equation (4.2-51) gives 


d 9b (t ' 



f) 




I 




dt ' 


f j- 

I(t - 

t* 

> ♦*[? 

i 


t-At ( 



\ 



+ w B/I (t) (t 

- 


dt 

• 

I- 1 

At 2 + 

1 

6 

[x 1 + “b/i 

<•>] 

-*At 

[i* 

At 

" ^eb^' 

t 

- At) 


(4.2-63) 



-At/r 

e 


Q 


b 


Substituting Equation (4.2-24) into this gives 

d Qb (t, t - At) * - J At A B/I (t)^ ^b/T ( t ) At 

+ A^(t, t - At) e" At/T ]o b 

Then, from Equation (4.2-56), 

t - At) * i A^jCt) At 

+ A*(t, t - At) e“ it/T j Q b 


(4.2-64) 


(4.2-65) 
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We divide the time interval from t' to t into n subintervals 
of length At = (t - t')/n as before. Then, from Equa- 
tions (2.1-8) and (2.2-13), with t^ ■ t* + iAt 


d(t i+1 , t') = 4>(t i + 1 , t.) d(t i# f) <t> (t i+1 , t ± ) + d(t i+1 , t L ) 

(4.2-66) 


or, using the partitioning of Equations (4.2-15) and (4.2-49) 


d ee (t i+i' 


t*> 


d 0b (t i+l' fc,) 


- ♦, 

00 (t i+l' fc : 

L> d 

CD 

CD 

ft 

H- 

+ 

^ee^i+i' 

V 

d 0b <fc i 

+ 

^©b^i+l' 

fc i> 

d 0b< fc i 

+ 

^eb^i+l' 

fc i> 

d bb (t i 

+ 

d ee (t i+i' 

t i ) 


" ^ee^i+i' t i 

L> ^ 

0b (t i' 

+ 

*eb (t i+l' 

V 

d bb (t i 

+ 

d 0b (t i+l' 

fc i> 



J ^ee (t i+i' t i ) 

t- > *eb (t i+l' fc i> 

fc, > *99 (t i+l' ‘i) 

t "> *9b (t i+l' *i> 

(4.2-67) 

* *bb (t i+l' t i ) 

fc *> ♦bb< t i +1 ' V 

(4.2-68) 
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Substituting Equations (4.2-19), (4.2-21), (4.2 24), 

(4.2-55), and (4.2-56) into these gives 


dA 


ee (t iti' t " > = a ee (t i' t ’ > + a eb (t i' t-) A b (t i+i' t i ) 


* R b (t Ul' V a 9b< t i' t ' ) 


(4.2-69) 


+ a; 


tp -wv f> A b< t i+1 , t ± ) 


* a ee^ t i+i* 


dA 


0b 


(t 1+1 , f) - d4b (t i' t ' ) * A b (t i + 1' V a bb<V t "> e 


-At/x 


(4.2-70) 


+ d 0b (t i+l' 


The matrices d£ Q (t, f) and d^(t, t‘) can be computed by re 
cursively applying Equations (4.2-69) and (4.2-70), starting 
from zero and using Equations (4.2—33), (4.2—57), (4.2 62), 

and (4.2-65) at each step. The matrices dg^t, t’) and 
d 0b(t # t’) can then be computed using Equations (4.2—55) 
and (4 . 2-56) . 

4.3 SENSOR MODELS 

The three— axis sensor modeled by ADEAS are 

• Gyroscopes 

• IR horizon sensor 

• Three-axis magnetometer 

• Fixed-head star tracker 

• Gimbaled star tracker 

• Digital Sun sensor 

• Analog Sun sensor 

• Gimbaled Sun sensor 
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For the three-axis stabilized spacecraft, ADEAS models the 
measurement function h(x, p') of Equation (2.1-13) in the 
following way: 


h(x, p’) = h(R S p*) (4.3-1) 

where R^ = unit vector, expressed in sensor coordinates, 

representing the known direction of an externally 
sensed object 

h = sensor-specific function of the unit vector R s 

The functional dependence of the measurement on the space- 
craft attitude and sensor alignment is through the unit ref- 
ac — 

erence vector R . Let A0 be the vector of small rotation 
angles representing the spacecraft attitude errors, and let 
<f> be the vector of small rotation angles representing the 
sensor alignment errors. The following equations assume 
that A© and 4> are expressed in radians. A0 is always a sub- 
set of the solved-for vector, whereas the analyst may des- 
ignate any component(s) of $ to be either solved for, 
considered, or ignored. The partial derivatives of the 
measurements with respect to A© and 4> are given as 


3v* = M_ Ml 

9A0 3R^ 3A0 


Ml _ M_ Ml , 4 o 

3<t> 3R b 3<t> 

AC 

3h/3R is sensor and model dependent and is described for 
each of the three-axis sensors in the subsequent sections; 
3R S /3A0 and 3R^/34> are derived below. 
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The unit reference vector in sensor coordinates can be ex 
pressed as 


R S - [As/b^W V - tA 


S/B 


] S B 


(4.3-4) 


where 


RS 


t A S/B^ 

[Ab/iI 

Rl 

RB 


the unit reference vector expressed in sensor 
coordinates 

the body— to— sensor rotation matrix 
the inertial— to— body rotation matrix 

the unit reference vector expressed in GCI 
coordinates 

the unit reference vector expressed in body 
frame coordinates 


The sensor alignment errors $ are errors in [A s/Q ] / and the 
spacecraft attitude errors A0 are errors in [A B/I ] . Thus, 
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considering only small errors in the attitude and sensor 

■'■VC 

alignments, the error in R can be expressed as 



The necessary partials are then 
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The partial derivatives 9h/9R^ in Equations (4.3-2) and 

A g 

(4.3-3) are evaluated assuming that the components of R can 
be independently varied. This is not strictly true because 
of the unit vector constraint |R^| = 1/ this constraint also 

A ^ - 

allows us alternate ways of expressing functions of R , for 
instance 



/s S 

Any resulting ambiguity in the partial derivatives 9h/9R is 
eliminated upon multiplication by 9R S /9A© or 9R S /9<£, so the 
partial derivatives 9y/9A0 and 9y/94> are unambiguous. The 
evaluation of 9R S /9A© is always performed because the atti- 
tude error vector is always a subset of the solved-for vec- 
tor, whereas 9R S /9<t> is computed only when one or more of the 
sensor alignment angles have been designated as solved for 
or considered. The formulation of h and its partial deriva- 
tives with respect to R^ is thus a key to the sensor— related 
error analysis computations. The expressions for h(R , p') 
and a complete list of all partial derivatives for each 
sensor are given below. 

For all sensors the sensor boresight is along the Z-axis of 
the sensor coordinate frame. 

4.3.1 GYROSCOPES 

Gyroscopes are modeled very differently from other sensors 
in ADEAS since gyroscope errors are regarded as part of the 
dynamic error model, as discussed in Section 4.2, rather 
than as measurement errors. Another major difference is 
that gyroscope data are assumed to be continuously available 
for attitude propagation, rather than being scheduled like 
other measurements. These differences can be summarized by 
saying that gyroscope data are used in a "model replacement 
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mode since they replace a full dynamic model of the space- 
craft motion incorporating environmental and control torques. 

Thus, in distinction to other sensor models, no measurement 
function h(x, p') will be derived for gyros. Rather, this 
section will provide a derivation of Equations (4.2-11) 
through (4.2-14). 

The raw gyro measurement w s is assumed to be related to 
the spacecraft angular velocity w B/I by 

Wg = M(w b/i - b - u 0 ) (4.3-7) 

where M is a matrix incorporating gyroscope alignments and 
scale factors, b is a vector of gyro biases, and Ug(t) is a 
white noise process. These quantities are not perfectly 
known, so the attitude determination is based on these meas- 
urements using a model 


to s = M*(a)* /]; - b*) (4.3-8) 

where M* and b* contain estimates of the alignments, scale 
factors, and biases. Thus, the estimate of the spacecraft 
angular velocity is 


“B/I = M * 1 M[c,) b/I " b " u 0 ] + b * (4.3-9) 

The matrix M* -1 M can be written as the sum of a symmetric 
and an antisymmetric matrix. 

M* -1 M = K + Ae (4.3-10) 
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where K is symmetric and Ae is skew-symmetric. It is assumed 
that M* is a good estimate of M; so K is close to the iden- 
tity matrix, and Ae is small. If K - I, the matrix I + Ae 
represents a small rotation given by a rotation vector Ae, 
the vector of gyro alignment errors. In the general case, 
we assume that Ae has this interpretation and that 

K = I + diag [ Ak] = I + diagtAk^, Ak^., Ak^] (4.3-11) 

where Ak is a vector of gyro scale factor errors. Nonzero 
off-diagonal elements of K would represent either higher- 
order effects in Ae, which are negligible, or shear-type 
misalignments of the gyro input axes (as opposed to the rigid 
misalignment of all input axes represented by Ae) , which we 
assume not to be present. 

Inserting Equations (4.3-10) and (4.3-11) int °_^ ua "_ T 
tion (4.3.9), and neglecting terms of order Akb , AkUg, 

Aeb, and AeUq gives 

= {I + diagtAk] + Ac) « B/I + b* - b - u Q (4.3-12) 
With Equation (4.2-9) written as 


Aw 


B/I 


= W 


B/I 


- w 


B/I 


(4.3-13) 


and with 


Ab = b* - b 


(4.3-14) 


4-35 


0450 


and 


Q = diag[co B/I ] 


(4.3-15) 


this can be written 


= Ab + Q Ak - Ae - Uq (4.3-16) 

which is identical to Equation (4.2-11). 

The remaining task is to specify the time dependence of the 
gyro biases, scale factors, and misalignments. The scale 
factor vector Ak and the misalignment vector Ac are assumed 
to be constant, as reflected in Equation (4.2-14). The bias, 
however, is assumed to obey the differential equation 


b = -b/r - u fa (4.3-17) 

where t is a correlation time and u b (t) is another white 
noise process. The processes Ug(t) and u^t) are assumed 
to be independent so that 


[ u e (t> 

^<f 

1 J 

II 

0 

CD 

cr> 

rf 

1 

f) 

(4.3-18) 

[v t} 

^(t*: 

)] = Q b «(t - 

c> 

(4.3-19) 

E 

c 1 

CD 

(T 

^(f)] = 0 


(4.3-20) 
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Equations (4.3-18) and (4.3-19) are identical to Equa- 
tions (4.2-53) and (4.2-54). The bias estimate obeys 

b* = -b*/T (4.3-21) 

with the same correlation time as b. Thus, Ab obeys the 
equation 

Ab = -Ab/x + u b (4.3-22) 

which is used in Equation (4.2-14). 

4.3.2 IR HORIZON SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The measurement model for the IR horizon sensor is shown in 
Figures 4-2 and 4-3, where the reference unit vector E is 
the Earth vector expressed in sensor coordinates. The out- 
puts of the IR horizon sensor are the Earth-in azimuth. 
Earth— out azimuth. Earth width, and Earth azimuth. Each 
azimuth is measured with respect to a user-defined, fixed- 
reference azimuth expressed in degrees from the X-axis of 
the sensor frame. The Earth vector in sensor coordinates is 
given by 

in n cos (<J> E + ^q) 

in n sin (<I> E + 4> 0 ) 
cos n 

where cos n = E S • Z S = E^ and 

0 since 0 < n < 180°. 
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Figure 4-3. Spherical Geometry for Earth IR Horizon Sensor 
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Thus the measurements are calculated as 


y* 


y* 


Earth width = h(E ) 


Q - 2 COS - 1 gQS p - COS n CPS Y 
I sin n sm y 


AC 

Earth azimuth = h(E ) = <I> 


(4.3-24) 


y* 

y* 


= tan 


e s 

-1 I - y 


- 4 >, 


Earth-in • h(E^) 


A IN " ^E 2 Q 


Earth-out = h(E S ) = A QUT = 4> E + \ Q 


(4.3-25) 

(4.3-26) 

(4.3-27) 


where E x , E y/ E z = x, y, and z components of Earth vector ex- 
pressed in the sensor frame# the super- 
script S being understood 

<l >0 = reference azimuth from X-axis of 
sensor frame (deg) 

p = Earth angular radius 
= sin -1 [(r + hj-)/R] 

r = Earth radius (km) 
h^ = IR tangent height (km) 

R = magnitude of the spacecraft position 
vector (km) 

y = sensor scan cone angle (deg) 

n = angle between sensor boresight and Earth 
vector (deg) 

The measurement model, in Equations (4.3-24) through 
(4.3-27), describes a conventional wheel-mounted horizon 
sensor. A very similar model describes the panoramic at- 
titude sensor (PAS) in the planar mode. Its unit reference 
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vector is also the Earth vector expressed in sensor 
coordinates. The differences are as follows: 

1. The sensor scan cone angle is 90 degrees. 

2. Only the Earth-in and Earth-out signals are sensed. 

3. The Earth-in and Earth-out azimuths are measured 
with respect to the X-axis of the sensor frame, 

i.e., the reference azimuth is 0 degrees. 

4. The sensor does not rotate continuously in one 
direction as does the sensor indexed scanner, but 
instead scans back and forth across a user-defined 
arc at a user-defined rate.- 

The direction of the scan is important for the computation 
of the measurements and their partial derivatives. When the 
sensor is scanning in the clockwise direction, the measure- 
ments are computed as 


= Earth-in 

AC 

= h(E S ) 

= A IN “ ^E + 2 n 

(4 

.3-28) 

= Earth-out 

AC 

= h(E S ) 

- A OUT ^E “ 2 Q 

(4 

.3-29) 


When the scanner is scanning counterclockwise. Equa- 
tions (4.3-25) and (4.3-26) are used for the Earth-in and 
Earth-out measurements, respectively. 

The measurement parameters that the user may select as either 
solved for or considered for the IR horizon sensor are 

1. The x sensor misalignment angle (deg) 

2. The y sensor misalignment angle (deg) 

3. The z sensor misalignment angle (deg) 

4. The Earth-in bias (deg) 

5. The Earth-in scale factor 
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6. The amplitude of the Earth-in periodic measurement 
error (deg) 

7. The Earth-out bias (deg) 

8. The Earth-out scale factor 

9. The amplitude of the Earth-out periodic measurement 
error (deg) 

10. The Earth width bias (deg) 

11. The Earth width scale factor 

12. The amplitude of the Earth width periodic measure- 
ment error (deg) 

13. The Earth azimuth bias (deg) 

14. The Earth azimuth scale factor 

15. The amplitude of the Earth azimuth periodic meas- 
urement error (deg) 

16. The sensor scan cone angle (deg) 

17. The fixed reference azimuth (deg) 

18. The IR tangent height (km) 

19. The distance from the spacecraft to Earth center 
(km) 

20. The Earth angular radius (deg) 

Note that even though the scan cone angle is always nomi- 
nally 90 degrees and the reference azimuth is always nom- 
inally 0 degrees for the PAS in the planar mode, both can 
have a specified uncertainty about their nominal values and 
hence can be designated as error parameters. 
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In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The Euler rotation sequence and Euler angles defin- 
ing the orientation of the sensor frame relative to 
the spacecraft body frame 

2. The frequency and phase angle of the Earth— in 
periodic measurement error 

3. The frequency and phase angle of the Earth-out 
periodic measurement error 

4. The frequency and phase angle of the Earth width 
periodic measurement error 

5. The frequency and phase angle of the Earth azimuth 
periodic measurement error 

6. The reference azimuth measured from the x-axis of 
the sensor frame (deg) 

7. The sensor scan cone angle (deg) 

8. The IR tangent height (km) 

9. The initial scan arc angle (deg) 

10. The time rate of change of the scan arc angle 

(deg/sec) 

11. The minimum scan arc angle (deg) 

12. The maximum scan arc angle (deg) 

Note that some of these parameters are not used for the PAS, 
and some are only used for the PAS. 

When scheduling the IR horizon sensor, the user may use one 
or two of the sensor outputs for the error analysis compu- 
tations. Since all outputs have units of degrees, the user- 
supplied value of the sensor white noise standard deviation 
for each measurement has units of degrees. 
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The errors in the Earth-in and Earth-out measurements are 
assumed to have equal standard deviations a A and to be un- 
correlated with one a no the i . Then, it follows from Equa- 
tions (4.3-26) through (4.3-29) that the errors in the 
Earth-width and Earth-azimuth measurements are uncorrelated 
and that their standard deviations are 



V 


fl 


(4.3-30) 


a Q 



(4.3-31) 


The correlations between the pairs (Q, A IN ) , (Q, a qut) ' 
(4> E/ A in ) , and (<I> E , a qut^ are i9 nore< i if any such pair is 
selected. 

Partial of the Earth azimuth wrt the unit Earth vector E: 



(4.3-32) 


(4.3-33) 


(4.3-34) 


Partial of the Earth azimuth wrt the sensor scan cone angle: 


9 ^ 

37 s ' 0 


(4.3-35) 
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Partial of the Earth azimuth wrt the Earth angular radius: 



(4.3-36) 


Partial of the Earth azimuth wrt the IR tangent height: 


J =0 (4.3-37) 

3h t 


Partial of the Earth azimuth wrt the spacecraft distance 
from Earth: 


3R 


0 


(4.3-38) 


Partial of the Earth azimuth wrt the bias b: 


3b 


1 


(4.3-39) 


Partial of the Earth azimuth wrt the scale factor k: 



h(E) 



(4.3-40) 


Partial of the Earth azimuth wrt the amplitude of the peri 
odic error a: 


3a 


sin (tot + ¥) 


(4.3-41) 


4-45 


0450 


Partial of the Earth width wrt the unit Earth vector E: 


ao. „ o 

8E x 


( 4 . 3 - 42 ) 


M. 

3E_ 


= 0 


( 4 . 3 - 43 ) 


_aa_ 


_ . 

3E 3(cos n) 

z 


= 2 


cos y - sin y cos 2 cot h 

. n 

sin y sin n sm — 


( 4 . 3 - 44 ) 


Partial of the Earth width wrt the sensor scan cone angle: 


9y " 2 


a 

cos y sin r| cos 2 - sin Y cos n 

. Q 

sm y sm n sm 


( 4 . 3 - 45 ) 


Partial of the Earth width wrt the Earth angular radius: 


3p 


= 2 


sin p 


Q 


sm y sin n sin — 


( 4 . 3 - 46 ) 


Partial of the Earth width wrt the IR tangent height: 


an _ in 

3h. " 3p 


r — 1 — 1 

I R cos pj 


( 4 . 3 - 47 ) 


Partial of the Earth width wrt the spacecraft distance from 
Earth : 


an _ _ in tan p 

8R “ dp R 


( 4 . 3 - 48 ) 
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Partial of the Earth width wrt the bias b: 


IS _ i (4.3-49) 

8b 


Partial of the Earth width wrt the scale factor k: 


Partial 
error a 


32 _ h ,s, . , cos -i f gaa p.- n 

3k ~ ^ cos I sin y sin n 


if the Earth width wrt the amplitude of the 


(4.3-50) 

periodic 


= sin (ut + Y) 
8a 


(4.3-51) 


In the following partial derivatives with a ± sign, the 
upper sign is to be used for the conventional IR scanner and 
for the PAS scanning in the counterclockwise direction, 
while the lower sign is used for the PAS scanning in the 
clockwise direction in accordance with Equations (4.3-28) 
and (4.3-29). 

Partial of the Earth-in angle wrt the unit Earth vector E: 


9A IN _ 9< *E 
3E x = 8E x 


(4.3-52) 


3A 

8E 


IN 

y 



(4.3-53) 


8A 

8E 


IN 

z 


= + 



(4.3-54) 
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Partial of the Earth-in angle wrt the sensor scan cone angle: 


dA 


IN 


3y 




(4.3-55) 


Partial of the Earth-in angle wrt the Earth angular radius: 


3A 


IN 


3p 


* ± (- 2 “) 


(4.3-56) 


Partial of the Earth-in angle wrt the IR tangent height: 


9a in _ . / 1 
3h fc ~ ± [- 2 8h t j 


(4.3-57) 


Partial of the Earth-in angle wrt the spacecraft distance 
from Earth: 


—IU _ . I 1 

3R ± l 2 3R ) 


(4.3-58) 


Partial of the Earth-in angle wrt the bias b: 


dA 


IE 


3b 


= 1 


(4.3-59) 


Partial of the Earth-in angle wrt the scale factor k: 


3A 


IN 


3k 


* h(i) = f E ± a) 


(4.3-60) 


Partial of the Earth-in angle wrt the amplitude of the 
periodic error a: 


3A 


IN 


3a 


= sin (tot + i|>) 


(4.3-61) 
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Partial of the Earth-out angle wrt the unit Earth vector E: 


3 A, 


OUT 


3E. 


3A 


OUT 


9E. 


9A OUT 
9E z =± 



(4.3-62) 

(4.3-63) 

(4.3-64) 


Partial of the Earth-out angle wrt the sensor scan cone 
angle : 


3A, 


OUT 


3Y 


= + 


(l 

\2 By) 


(4.3-65) 


Partial of the Earth-out angle wrt the Earth angular radius: 


3 A, 


OUT 


3p 




(4.3-66) 


Partial of the Earth-out angle wrt the IR tangent height: 


9A OUT (1 90 - \ 
3h t = ± [2 3h fc ) 


(4.3-67) 


Partial of the Earth-out angle wrt the spacecraft distance 
from Earth: 


9A OUT /l 90 \ 

9R ~ ± \ 2 9R y / 


(4.3-68) 
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Partial of the Earth-out angle wrt the bias b: 


8a out _ 
9b 


(4.3-69) 


Partial of the Earth-out angle wrt the scale factor k: 

h(E> - * E ± (A Q) (4.3-70) 

Partial of the Earth-out angle wrt the amplitude of the 
periodic error a: 


9A OUT 

8a 


= sin (cot + y) 


(4.3-71) 


4.3.3 THREE -AXIS MAGNETOMETER MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The three-axis magnetometer computes the three independent 
components of the Earth's magnetic field at any time and 
position, using the geomagnetic reference field updated to 
epoch 1985.0. The measurement is the magnetic field vec- 
tor M, which can be written M * MU, where M is the magnitude 
of the magnetic field, and U is a unit vector in the direc- 
tion of the field. The external reference vector is the 
unit vector U: 


y* = x component 


y* = y component 


y* = z component 


h(U) = = MU^ 

h(U) = = MUy 

h(U) = M S = MU S 

z z 


(4.3-72) 


(4.3-73) 


(4.3-74) 
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The measurement parameters, which the user may select as 
either solved for or considered, for the three-axis magnetom- 
eter are 


1 . 

The x sensor misalignment angle (deg) 


2. 

The y sensor misalignment angle (deg) 


3 . 

The z sensor misalignment angle (deg) 


4 . 

The x component bias (Teslas) 


5. 

The x component scale factor 


6. 

The amplitude of the x component periodic 
ment error (Teslas) 

measure- 

7. 

The y component bias (Teslas) 


8. 

The y component scale factor 


9 . 

The amplitude of the y component periodic 
ment error (Teslas) 

measure- 

id. 

The z component bias (Teslas) 


11. 

The z component scale factor 


12. 

The amplitude of the z component periodic 
ment error (Teslas) 

measure- 


In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The Euler rotation sequence and Euler angles defin- 
ing the orientation of the sensor frame relative to 
the spacecraft body frame 

2. The frequency and phase angle of the x-component 
periodic measurement error 

3. The frequency and phase angle of the y-component 
periodic measurement error 
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4. The frequency and phase angle of the z-component 
periodic measurement error 

When scheduling the three-axis magnetometer, the analyst may 
use any or all of the sensor outputs for the error analysis 
computations. Since all three outputs are components of a 
magnetic field vector, the sensor white noise standard 
deviation for each measurement has units of Teslas. 


Partial of the x component wrt the unit magnetic field vec- 
tor U: 



(4.3-75) 



(4.3-76) 



(4.3-77) 


Partial of the x component wrt the bias b: 



(4.3-78) 


Partial of the x component wrt the scale factor k: 


8M 
£ 

8k 


~ h(U) 



(4.3-79) 
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Partial of the x component wrt the amplitude of the periodic 
error a: 


9M 

- — = sin (cot + Y) 


( 4 . 3 - 80 ) 


Partial of the y component wrt the unit magnetic field vec 
tor U: 




( 4 . 3 - 81 ) 


( 4 . 3 - 82 ) 


( 4 . 3 - 83 ) 


Partial of the y component wrt the bias b: 



( 4 . 3 - 84 ) 


Partial of the y component wrt the scale factor k: 


Partial of the y 
error a: 



component 


« h(U) = M y 

wrt the amplitude 


of 


( 4 . 3 - 85 ) 
the periodic 



sin (cot + t|r) 


( 4 . 3 - 86 ) 
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Partial of the z component wrt the unit magnetic field vec- 
tor “0: 





M 


Partial of the z component wrt the bias b: 



(4.3-87) 


(4.3-88) 


(4.3-89) 


(4.3-90) 


Partial of the z component wrt the scale factor k: 


3M 

* h(U) = M z (4.3-91) 


Partial of the z component wrt the amplitude of the periodic 
error a: 


3M 

= sin ( wt + i|i) (4.3-92) 

4.3.4 FIXED-HEAD STAR TRACKER MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The fixed-head star tracker projects the x and y coordinates 
of the reference unit star vector S S onto the U-V plane, as 


shown in Figure 4-4. The measurements are computed as fol- 
lows : 


y« . 0 - h(S S ) = - (4.3-93) 

y* * V = h(S S ) - - (Sy/S^) (4.3-94) 

ADEAS allows the fixed-head star tracker to have either a 
conical or pyramidal field of view. The star is visible to 
the conical sensor only if the angle between the boresight 
and S 8 is less than the user-specified conical half angle. 
Similarly, the star is visible to the pyramidal star tracker 

only if 

Pyramidal x half angle > 

and 



. -l /„s/ c ,s\ 

Pyramidal y half angle > 

tan [ S Y / S z ) 



The measurement parameters, which the user may select as 
either solved for or considered, for the fixed— head star 
tracker are 

1. The x sensor misalignment angle (deg) 

2. The y sensor misalignment angle (deg) 

3. The z sensor misalignment angle (deg) 

4. The U bias 

5. The U scale factor 

6. The amplitude of the U periodic measurement error 


•W> 


\ 
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7. The V bias 

8. The V scale factor 

9. The amplitude of the V periodic measurement error 

In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The Euler rotation sequence and Euler angles defin- 
ing the orientation of the sensor frame relative to 
the spacecraft body frame 

2. The frequency and phase angle of the U periodic 
measurement error 

3. The frequency and phase angle of the V periodic 
measurement error 

4. The designation of the sensor field of view as 
either conical or pyramidal 

5. The conical half angle or the pyramidal x and y 
half angles 

When scheduling the fixed-head star tracker, the analyst may 
use either or both of the sensor outputs for the error anal- 
ysis computations. Since both outputs are dimensionless 
projections of a unit vector, the sensor white noise 
standard deviation has units of radians. 

A 

Partial of the U component wrt the unit star vector S: 



(4.3-95) 
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0 


(4.3-96) 


ILL _ 
3Sy ' 


2U_ 

3S_ 


L 

s 2 . 


Partial of the U component wrt the bias b: 


3S1 _ i 
3b ~ 1 


Partial of the U component wrt the scale factor k: 


5IZ w$\ _ 

3k ~ h(S) " S 


Partial of the U component wrt the amplitude of the 
error a: 


= sin (wt + »|i) 


Partial of the V component wrt the unit star vector 


2V_ 

3S 


» 0 


3v_ 

as 


i_ 

s 


$v_ 

as 


5* 
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(4.3-97) 

(4.3-98) 

(4.3-99) 

periodic 

(4.3-100) 

/V 

S: 

(4.3-101) 

(4.3-102) 

(4.3-103) 
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Partial of the V component wrt the bias b: 


5V _ 1 <4,3-104) 

9b " 

Partial of the V component wrt the scale factor k: 

$**h<S)-^ (4.3-105) 

z 

Partial of the V component wrt the amplitude of the periodic 
error a: 


22 

8a 


sin (oat + 


(4.3-106) 


4.3.5' GIMBALED STAR TRACKER MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 


The gimbaled star tracker uses azimuth (4>) and elevation (9) 
gimbals to point a focal plane at a star. Thus, if the star 
unit vector is S^ in the sensor frame, it is % in the focal 
plane frame, where (see Figure 4-5) 



"sin 

© 

0 

-cos 0 ” 


" COS 

4 > 

sin 4> ( 

3 1 

a FP = 

0 


1 

0 


-sin 

4 > 

cos <j> 0 


cos 

0 

0 

sin © _ 


0 


0 









(4 


sin 

0 

(cos 

4> S + sin 
z 

♦ S y ) 

- 

COS 0 S z 



.cos <t> S - sin 4> S 

y a 

cos 0 (cos <j> S + sin 4> S ) + sin 0 S z 

JL Jr 


The superscript S has been omitted on the right side for 
convenience. A superscript S should be assumed where no 
superscript appears. 
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The gimbals are nominally controlled to center the star in the 
focal plane, i.e., to make S FP = [0 0 1] , and the gimbal 

angles are the measurements. Thus, the measurements are 
computed as 

y* = azimuth E <t> = h(S^) = tan ^ (4.3-108) 

Y* = elevation E 0 = h(S^) = sin (^z) (4.3-109) 

The azimuth equation is to be interpreted as 


(4.3-110) 


(4.3-111) 

The field of view of the gimbaled star tracker is specified 
by a minimum and maximum elevation and a range of azimuths. 
The azimuth range is specified by its midpoint and full 
width; if a width of more than 360 degrees is specified, the 
field of view includes all azimuth values. 

The sensor noise modeling must allow for imperfect centering 
of the star image in the sensor focal plane. Thus, if we let 

4> = tan -1 (S /S ) + A4> (4.3-112) 

y x 

0 = sin ^ (S ) + A0 (4.3-113) 

z 



• - '• \ 
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(4.3-114) 


£ FP 


[1 - e_ - e' 


we find, to first order in the small quantities A<£, A0, 
e x , and e y 


A4> = -e 




(4.3-115) 


A0 = e x 


(4.3-116) 



To these gimbal errors due to focal plane centering errors, 

we must add intrinsic errors in the gimbal readouts. We 

take the mean squares of A<(> and A0, assuming equal variances 
2 

o F p in the x and y focal plane errors and assuming the 
different errors to be uncorrelated. This gives 



variance in = 



(4.3-117) 


2 2 

variance in © = Oq + o Fp (4.3-118) 

2 2 

where and Oq are the variances of the gimbal readout 
errors and o^, Oq, and a Fp are all given in degrees, the 
units of the measurements. 

The measurement parameters, which the user may select as 
either solved for or considered, for the gimbaled star 
tracker are 

1. The x sensor misalignment angle (deg) 

2. The y sensor misalignment angle (deg) 
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3. The z sensor misalignment angle (deg) 

4. The azimuth bias (deg) 

5. The azimuth scale factor 

6. The amplitude of the azimuth periodic measurement 
error (deg) 

7. The elevation bias (deg) 

8. The elevation scale factor 

9. The amplitude of the elevation periodic measurement 
error (deg) 

In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The Euler rotation sequence and Euler angles defin- 
ing the orientation of the sensor frame relative to 
the spacecraft body frame 

2. The frequency and phase angle of the azimuth pe- 
riodic measurement error 

3. The frequency and phase angle of the elevation pe- 
riodic measurement error 

4. The minimum and maximum azimuth and elevation de- 
fining the field of view 

When scheduling the gimbaled star tracker, the analyst may 
use either or both of the sensor outputs for the error 
analysis computations. 

Partial of the star azimuth 4> wrt the unit star vector S: 

(4.3-119) 

2 r,2 


94 > 

9S x = " 
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2 2 

S x + S y 


(4.3-120) 


Partial of the star azimuth wrt the bias b: 


(4.3-121) 


M _ 

9b “ 1 


(4.3-122) 


Partial of the star azimuth wrt the scale factor k: 


3k* h <*> = tan_1 S 


( 4 . 3/*-123 ) 


Partial of the star azimuth wrt the amplitude of the 
periodic error a: 


a* 

= sin (cut + \|») 


(4.3-124) 


Partial of the star elevation © wrt the unit star vector S: 


(4.3-125) 


(4.3-126) 


1 

cos 9 




(4.3-127) 


4-64 


0450 



Partial of the star e 


levation wrt the bias b: 


IS _ i (4.3-128) 

8b 


Partial of the star elevation wrt the scale factor k: 


« h(S) = sin 1 (S y ) ( 4 - 3 129 ) 


Partial of the star elevation wrt the amplitude of the per 
iodic error a: 


■§® = sin (ut + (4.3 130) 

3a 

If S » 1, the quantity + S y should be given a small, 
finite positive value to avoid division by zero. 

4.3.6 DIGITAL SUN SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The digital Sun sensor provides measurements of the angle a 
between the sensor boresight and the projection of the unit 
Sun vector S S in the Y-Z plane and the angle (3 between the 
sensor boresight and the projection of the unit Sun vector 
in the X-Z plane, as shown in Figure 4-6. The sensor out- 
puts N and N^ are related to these angles by Equa- 
tions (7-38) and (7-39) of Reference 3: 

a = tan -1 <S y /S z ) = A g + tan -1 [A x + AgNJ + A 3 sin (A 4 N a 

(4.3-131) 

+ A 5 ) + A 6 sin (A ? N a + A q )] 

(3 = tan -1 <S X /S z ) = B g + tan -1 [B 1 + + Bg sin (B 4 N R 

(4.3-132) 

+ B 5 ) B 6 sin (B 7 Ng + Bg)] 
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where A 


provided’ by the manufacturer, and where the superscript S is 
to be understood on the components of the Sun unit vector 
where no superscript is shown. 


Ag and B 1 , 


Bg are calibration constants 


The constants Aj, A 6> B 3 , and B 6 are very small. If A 3 

k c are neglected, Equation (4.3-131) gives 
6 


and 


A 2 N a “ "a 


- tan (a - A q ) - A. 


S - S„ tan A- 
y z 2 . 

s y tan Ag + S z 


- A. 


(4.3-133) 


This expression for can be substituted into the sinusoidal 
terms on the right side of Equation (4.3-131), and a similar 
procedure can be used in Equation (4.3-132) to obtain meas- 
urement equations valid to first order in A 3 , A 6 , B 3 , and B 6 . 
Thus, the measurements are computed as 

y* = N a = h(s , A^, * • • • ) 

= [n a " A 3 sin (A 5 + A 4 n a /A 2 > 


(4.3-134) 

- A 6 sin (A g + A 7 n a /A 2 )]/A 2 


y* = Njj = h( , B^, • • • • ®g) 

= [n 0 - B 3 sin (B 5 + B 4 n fJ /B 2 ) - B g sin (B Q 
where n is given by Equation (4.3-133) and 

CL 

= tan (B - B g ) - B 1 

S - tan B q 
S x tan B 9 + S z 1 


(4.3-135) 

+ B ? n 0 /B 2 )]/B 2 


(4.3-136) 
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In ADEAS, it is assumed that the Sun sensor outputs are in 
units of degrees. 

ADEAS allows the Sun sensor to have either a conical or a 
pyramidal field of view. The Sun is visible to the conical 
sensor only xf the angle between the boresight and S is 
less than the user-specified conical half angle. Similarly, 
the Sun is visible to the pyramidal Sun sensor only if 


and 


Pyramidal x half angle > 


tan 



Pyramidal y half angle > 


tan 



For either configuration, ADEAS also determines whether the 
sensor is occulted by the Earth. 


The measurement parameters, which the user may select as 
either solved for or considered, for the two-axis digital 
Sun sensor are 


1 . 
2 . 
3. 
4 . 


5. 


6 . 


7. 


The x sensor misalignment angle (deg) 

The y sensor misalignment angle (deg) 

The z sensor misalignment angle (deg) 

The N a calibration constants A^ , A 3 , A 4 , A g/ and A ? 
(dimensionless) 

The N a calibration constant A 2 (deg -1 ), nominal 
value ir/180. 

The N calibration constants A c , A Q , and A- (deg) 
cl o o y 

The N bias (deg) 
a 
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8. 

The N scale factor 
oc 



9. 

The N amplitude of the periodic 

OC 

measurement 

error 


(deg) 



10. 

The Np calibration constants B 1 , 

B 3' B 4 ' B 6 ' 

and B 7 


(dimensionless) 



11. 

The N q calibration constant B 2 (deg 1 ), nominal 


value ir/180 . 



12. 

The Np calibration constants B,., 

B q , and Bg 

(deg) 

13 . 

The bias (deg) 



14. 

The Njj scale factor 



15. 

The N £ amplitude of the periodic 

measurement 

error 


(deg) 




In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The Euler rotation sequence and Euler angles defin- 
ing the orientation of the sensor frame relative to 
the spacecraft body frame 

2. The frequency and phase angle of the a periodic 
measurement error 

3. The frequency and phase angle of the 3 periodic 
measurement error 

4. The designation of the sensor field of view as 
either conical or pyramidal 

5. The conical half angle or pyramidal x and y half 
angles 

When scheduling the two-axis Sun sensor, the user may use 
either or both of the sensor outputs for the error analysis 
computations. Since the outputs have units of degrees, the 
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value of the sensor white noise standard deviation for 
measurement has units of degrees. 

A 

Partials of N wrt the unit Sun vector S: 

CL 

aN ct /9S x " 0 (4 


3N /8S 

(X 


y 


3N /3S 
a 


where 


z 


(3N a /3n a ) s z (l + tan 2 A g )/(S z + S y tan A g ) 2 

(4 

(3N a /3n a ) S y (l + tan 2 A g )/(S z + S y tan A g ) 2 

(4 


3N /3n = [A_ - A,A. cos (A_ + A. n /A.,) 

a a 2 34 5 4a 2 

*" cos (Ag + Ay ) 1 /A^ 

Partials of N wrt the calibration constants: 
a 

3 V 3A 1 - - 3 V 3n a 

3 V 3A 2 = -N a /A 2 + (n a /A^)<l - A 2 3N a /3n a ) 

3N /3A- = -A” 1 sin (A- + A. n /A-) 
a 3 2 5 4 a 2 

3N a /8A 4 = -( A 3 n a /A 2 ) cos (A 5 + A 4 n a /A 2> 
SN^/SAg = -(A 2 /A 2 ) cos (A 5 + A 4 n a /A 2 ) 


(4. 

(4. 

(4. 

(4. 

(4. 

<4. 
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,3-138) 

3-139) 

3-140) 

3-141) 

3-142) 

3-143) 

3-144) 

3-145) 
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3N 


/3A-. = -A- 1 sin (Ag + A ? n a /A 2 ) (4.3 146) 


a 6 2 


3N /3A_ = 
oc / 


"( A 6 n a /A 2 ) COS (A 8 + A 7 n a /A 2 } (4 * 3 14?) 


3N 


/3A fl = -(Ag/A 2 ) COS (Ag + A ? n a /A 2 ) (4.3 148) 


a 8 


3N a /3A 9 = -OV 3 n a)( S y + s z)( 1 * ta " 2 A 9y <S * + s y 1 


an A g ) 
4.3-149) 


Partial of N wrt the bias b s 

OC ^ 


3N /3b = 1 

a a 


(4.3-150) 


Partial of N wrt the scale factor k : 

OC w 


3N a /3k a « h|s S , A 1# ..., A g j = (4.3 151) 


Partial of N wrt the amplitude of the periodic error: 

OC 


3N /3a = sin (ait + Y) 

a a 


(4.3-152) 


Psrtials of N|j wrt the unit Sun vector S. 


9N (3 /3S x = < 9N 3 /8n |3> S z( 1 + tan B 9 


* 


S x tan B g ) 


(4.3-153) 


8 V 9S y = 0 


(4.3-154) 
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3N b /3S z = -<3N R /3n 3 ) S^l + tan 2 B^ /(S z + S % tan B g ) 


( 4 . 3 - 155 ) 


where 


3N B /3n B - t®2 - B 3 B 4 cos < B 5 + B 4 V B 2> 

- B 6 B 7 COS (B g 4 B ? n B /B 2 )]/B a 
Partials of N R wrt the calibration constants: 

= -3N (3 /3n R 

a V 3B 2 ’ -V B 2 + (V B 2) (1 - B 2 aN B /an B ) 
3 V 3B 3 = -B- 1 sin (B 5 + B 4 n B /B 2 ) 

a V 3B 4 * ~( B 3 b B /B 1) C0S < B 5 * B 4 V B 2> 

3N r /3B 5 = -(B 3 /B 2 ) cos (B 5 + B 4 n R /B 2 ) 

3 V 3B 6 * -*? sin < B 8 - B 7 V B 2> 
3 V 3B 7 - -( B 6 "B /B 2) C0s < B 8 + B 7 V B 2> 

3N r / 3B 8 = -(Bg/B 2 ) cos (Bg + B ? n^/B^ 

3 V 3B 9 = -(3M B /3n B )(s a 4 S 2 )(l 4 tan 2 * 9 )I<S Z 4 S;c 


( 4 . 3 - 156 ) 

( 4 . 3 - 157 ) 

( 4 . 3 - 158 ) 

( 4 . 3 - 159 ) 

( 4 . 3 - 160 ) 

( 4 . 3 - 161 ) 

( 4 . 3 - 162 ) 

( 4 . 3 - 163 ) 

( 4 . 3 - 164 ) 

tan B g ) 2 
( 4 . 3 - 165 ) 
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Partial of N R wrt the bias b R : 


aNjj/Sb 



(4.3-166) 


Partial of N^ wrt the scale factor k^: 


9 V 3k B 


« h(S S , B 1 .B 9 ) = Hg (4.3-167) 


Partial of N„ wrt the amplitude of the periodic error 


aN^/aa^ = sin (u>t + ¥) . 


(4.3-168) 


4.3.7 ANALOG SUN SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The geometry of the one-amis analog Sun sensor is shown in 
Figure 4-7. The z-amis is chosen to be the boresight of the 
sensor. The unit vector S S is the direction t 9 the Sun in 
sensor coordinates. The one-amis CSS has two simple photo- 
cell "eyes" whose outputs, f 1 and f 2 , are proportional to 
the cosine of the angle between the spacecraft to Sun vector 
and the normal to the eye surface. The sensor output f css 
is zero unless both eyes are illuminated, in which case the 
output is the difference of the two eye outputs. Thus the 
measurement is computed as: 


y* = h(S^) 

" f CSs( S x' 



(4.3-169) 

= l' 1 " 

f 2 if S S . 

/\ 

u l' 

* C 1 

AC 

and U 2 > C 2 

I ° 

if § s • 

/V 

u l' 

A 

O 

or S & • U 2 < C 2 
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where 


«! * 1 a §S * Ui - — s 


l|s® cos a + sin. aj (4.3-170) 


1 + b 
f 2 2 S 


- 1 1 b 2 / r S 


2- S S • U 2 = ~ 2 s ^ ( S z COS a " S x Sin a ) < 4 - 3 " 171) 


( £ 


The parameters b x and b 2 are scale factor errors; C 1 and C 2 
are the cosines of user-supplied fields of view of the two 
eyes, and U 1 and U 2 are unit vectors specifying the normals 

of the eye surfaces; 


A 



sin a 
0 

cos a 


A 



-sin a 
0 

cos a 


(4 . 3-172a) 


(4 . 3-172b) 


The calibration factor s is nominally equal to sin a. 

Thus, with the Sun vector in sensor coordinates written as 


S S 


cos © sin $ 
sin 9 
cos 0 cos 


(4.3-173) 


and with b^ = b 2 
sensor output is 


AC 

0, s 5 


• u, > c 


1 ' 


and S S 


* U 2 * C 2 ' 


the 


£ CSS - cos 6 sin + 


(4.3-174) 
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This output is nominally equal to $ in radians if both <J> 
and 9 are small angles, which is if the Sun is close to 
the sensor boresight. Note that errors in s do not have to 
be considered independently, since they are equivalent to 
errors in the scale factor k. 

The measurement parameters, which the user may select as 
either solved for or considered, for the analog Sun sensor 
are 

1. The measurement bias (radians) 

2. The measurement scale factor 

3. The amplitude of the measurement periodic error 
( radians) 

4. The half-angle between the eye normals, a (deg) 

5. The individual eye scale factor errors, b.^ and b 2 

In addition to specifying a name and uncertainty for each of 
the measurement parameters designated as either solved for 
or considered, the analyst must also provide the following: 

1. The frequency and phase factor of the periodic 
measurement error 

2. The half-angle between the eye normals, a (deg) 

3. The individual eye scale factor errors, b.^ and b 2 

(nominally zero) 

4. The cosines of the eye fields of view, and C 2 

When scheduling the analog Sun sensor, the analyst does not 
have a choice of measurements to use because there is only 
one . 

n computing the partial derivatives, s is set equal to 
sin a after carrying out the differentiations. The values 
given by Equations (4.3-175) through (4.3-183) are for the 
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case that both eyes are illuminated. If either 

§ s . 0 < C x or S S • U 2 < C 2 , the partial derivatives are all 

identically zero. 

Partial wrt the Sun vector S (superscript S to be under- 
stood) : 


3f CSS 2 + b l + b -2 
9S_ “ 2 


(4.3-175) 


9£ css = 0 

3S y - 


(4.3-176) 


ifcss . b l 'Ji 

8S_ " 2 


cot a 


(4.3-177) 


Partial wrt the half angle a: 


3f ___ 2 + b, + b 9 

CSS 1 2. 

8a 2 


S cot a 

x 



(4.3-178) 


Partials wrt the scale factors b^ and b 2 


^ - * (S, cot a ♦ 8,, 
^icss . . 1 (S z cot a - s x ) 


Partial wrt the bias b: 


(4.3-179) 


(4.3-180) 


a£ CSS _ •, (4.3-181) 

8b 
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Partial wrt the scale factor k: 


3f 


CSS 


dk 


= f 


CSS 


(4.3-182) 


Partial wrt the amplitude of the periodic error 


3f 


CSS 


da 


= sin (wt + <(>) 


(4.3-183) 


4.3.8 GIMBALED SUN SENSOR MEASUREMENT MODEL AND PARTIAL 
DERIVATIVES 

The gimbaled Sun sensor measurement model and partial deriv- 
atives are identical to those for the gimbaled star tracker, 
described in Section 4.3.5, except that the Sun unit vector 
is used in place of the star unit vector. The numerical 
values of parameters in the models may differ, of course. 
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SECTION 5 - REFERENCE SYSTE MS AND VECTORS 


5 . 1 COORDINATE SYSTEMS 

5.1.1 GEOCENTRIC INERTIAL (GCI) COORDINATE SYSTEM 

The GCI coordinate system, illustrated in Figure 5-1, defines 
a fixed set of axes in inertial space. The Zj axis^is paral- 
lel to the Earth's spin vector or north pole. The Xj axis 
is along the vernal equinox, the Earth to Sun direction when 
the Sun crosses the equatorial plane. The Yj axis is given 
by the vector cross product x X^. . In summary, 


/\ 



vernal equinox 
Z 1 x Xj/IZj x Xjl 
Earth spin vector 


(5.1-1) 


Z GCI 



Figure 5-1. 


Definition of the GCI Coordinate System 
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The GCI coordinate system can be true-of-date, mean-of-date, 
or mean-of-epoch for some fixed epoch. Only mean-of-epoch 
coordinate systems are truly inertial, but the effects of 
the motion of true-of-date and mean-of-date coordinate sys- 
tems are negligible for attitude error analysis computations. 
For high accuracy applications, all reference vectors (stars. 
Sun, Moon, spacecraft) should be expressed in the same iner- 
tial coordinate system. 

5.1.2 EARTH-CENTER POINTING COORDINATE SYSTEM 

/\ /N y\ 

The Earth-center pointing coordinate system (X_., Y„, Z„) is 

r* it E 

defined by the spacecraft orbital position and velocity. 

These give two orthogonal axes: E ="e/|e|, where f is the 

spacecraf t-to-Earth vector; E x N, where N is the orbit nor- 
mal vector. The vector +E can be chosen as one coordinate 
axis of the Earth-center pointing system, and +E x N as a 
second axis. Then the third axis is defined to give a right- 

/v /s /N 

hand orthogonal triad X„, Y_, Z„. 

£j E* XL 

The transformation matrix from GCI coordinates to Earth- 
center pointing coordinates is 





a r/i 



(5.1-2) 


^I -^1 /'•T 

where X E , Y E , and Z £ are 3x1 column vectors containing the 



Z £ , respectively. 
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5.1.3 LOCAL VERTICAL POINTING COORDINATE SYSTEM 

/\ ^ /\ . 

The local vertical pointing coordinate system (X^, Y^, Z L ) 
is defined in terms of the local vertical pointing ($ LV ) and 
orbit normal (N) unit vectors. The local vertical pointing 
unit vector components in GCI coordinates are given by 


P 1 

*LV 


-cos cos X 
-cos <j> sin x 
-sin $ 


(5.1-3) 


where x = geodetic longitude 


= arc tan [R y /R x ] 


4> = geodetic latitude 


= arc tan \ R, 


(1 - e ) 


4 


r2 + R J 

x y 


where e = Earth eccentricity and R = -E 1 is the Earth-to 
spacecraft vector resolved in GCI coordinates. 

A A A A 

The user specifies which unit vector (-X^, + x i/ -Y L' +Y L' 

-Z , +Z ) is defined as the pointing axis ( R LV > and which 
unit vector is defined as the axis P LV * N. The third axis 
is defined to give a right-hand orthogonal triad X L , Y^, 

The transformation matrix from GCI coordinates to local ver- 
tical pointing coordinates is 


-1 


(5.1-4) 


(5.1-5) 
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/NJ /SJ /Nj 

where X L , Y^, and are 3x1 column vectors continuing the 
components along the GCI axes of the vectors Y^, and Z L , 

respectively. 

5.1.4 SUN-POINTING COORDINATE SYSTEM 

A A A 

The Sun-pointmg coordinate system (X g/ Y g/ Z g ) is defined 
in terms of the following unit vectors: 


S = Spacecraf t-to-Sun Unit Vector 


R ES R l 

r es ~ R I 


XN /-V 


U = S x N t 


where R gg = Earth-to-Sun vector 

Rj = Earth-to-spacecraf t vector (= -E) 

N £ E GCI z-axis (North Pole) 

/\ 

The vector +S is chosen as one coordinate axis of the Sun- 

pointing system, and +U as a second axis. Then the third 

axis is defined to give a right-hand orthogonal triad X„, 

/N /N S 

V V 

The transformation matrix from GCI coordinates to Sun- 
pointing coordinates is 


(n r 


(5.1-9) 
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where x|, Y*, and z* are 3 x 1 column vectors^containing the 
components along the GCI axes of the vectors X g , Y g , Z g , 
respectively . 

5.2 SPACECRAFT EPHEME RTS GENERATION 

ADEAS has two options for spacecraft epheraeris generation: 
internal ephemeris generation by numerical integration and 
reading of an externally-generated ephemeris file. 

5.2.1 NUMERICAL INTEGRATION 

The Cartesian position r and velocity v of the spacecraft in 
an inertial coordinate system, obey the coupled ordinary 
differential equations 


dr/dt = v (5.2-1) 

dv/dt = a (5.2-2) 

where a is the acceleration of the spacecraft. The next 
three sections specify the algorithms for computing a. 

Then, the last two sections present the numerical integra- 
tion and interpolation algorithms. 

5. 2. 1.1 Earth Potential Including J 2 Term 

The force per unit mass on the spacecraft due to the Earth s 
gravity is 


= -px/r 3 {l - 

(3/2) 

J 2 (R/r ) 2 

[5(z/r) 2 

- 1]} 

(5.2-3) 

= -py/r 3 (l - 

(3/2) 

J £ (R/r ) 2 

[ 5 (z/r) 2 

- 1]} 

(5.2-4) 

= -]iz/r 3 { 1 - 

(3/2) 

J 2 (R/r) 2 

[5(z/r) 2 

- 3]} 

(5.2-5) 
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where 


u = 3.986005 x 10 5 km 3 /s 2 (5.2-6) 

is the default gravitational constant of the Earth 

J 2 = 1.08263 x 10 -3 (5.2-7) 

is the default geopotential coefficient for the Earth's 
oblateness , 


R = 6378.140 km (5.2-8) 

is the default equatorial radius of the Earth, and x, y, and 
z are the GCI components of the spacecraft position vector. 

Higher order effects of the nonsphericity of the Earth are 
not modeled in ADEAS. If these terms are important for an 
application, an externally-generated spacecraft ephemeris 
must be used. 

On user option, the terms can be omitted. The central 

Z — 3 

body term of the Earth, -pr/r , is always included. 

5.2. 1.2 Atmospheric Drag Effects 

Atmospheric drag acceleration is modeled as a drag force in 
the direction of the relative wind vector acting on a satel- 
lite of constant surface area. The velocity of the satel- 
lite relative to the atmosphere is computed in the inertial 
coordinate system by subtracting the motion of the atmos- 
phere, assumed to rotate with the Earth, from that of the 
satellite : 


V 


rel 


= r - to x r 


(5.2-9) 
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The Earth’s rotation vector, w, is directed along the 
Earth's instantaneous spin axis with a magnitude equal to 
the rotation rate of the Earth and components (oi r w 2 , oj 3 ) . 

Ignoring the effects of precession and nutation, the Z-axis 
is aligned with the north polar spin axis such that and 

^ -i t e o c a o a v in rafl/SPC 


0) 


a iiyntiu vuui * * 

ace equal to zero, and a> 3 = 7.29211585494 x 10* red/aec 


is constant. Therefore, in this approximation, Equa- 
tion (4-12) reduces to 


rel 


x + co 3 y 
y - co 3 x 
z 


( 5 . 2 - 10 ) 


For the case of a spherical satellite, the atmospheric drag 
acceleration is computed as 


= - 


1 

2 


m 


P 0 <»> v rel |V rel' 


(5.2-11) 


where C D = aerodynamic force coefficient, which is an ad 
justable parameter 

A = cross-sectional area of the satellite 
m = mass of the satellite 

p (h) = altitude density function computed from the 
0 atmospheric drag model 

Nominally, for a spherical satellite, the aerodynamic force 
coefficient, C D , is equal to 2.0. In ADEAS, the altitude 
density function, p Q (h), is modeled using a Harns-Pnester 
atmospheric model. Harris and Priester determined the phys- 
ical properties of the upper atmosphere theoretically by 
solving the heat conduction equation under quasi-hydrostatic 
conditions. Approximations for fluxes from the extreme 
ultraviolet and corpuscular heat sources were included, but 
the model averages the semiannual and seasonal-latitudinal 
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variations and does not attempt to account for the extreme 
ultraviolet 27-day effect. 


The atmospheric model presented here is a modification of 
the Harris-Priester concept. The modification attempts to 
account for the diurnal bulge by including a cosine varia- 
tion between a maximum density profile at the apex of the 
diurnal bulge (which is located approximately 30 degrees 
east of the subsolar point) and a minimum density profile at 
the antapex of the diurnal bulge. Discrete values of the 
maximum- and minimum-density altitude profiles, correspond- 
ing to mean solar activity, are stored in tabular form as 

^max^i^ an< ^ Pmin^i^' res P ec tively. Different maximum and 

minimum profiles are available for different levels of solar 

activity. Exponential interpolation is used between entries; 

i.e., the minimum and maximum densities, p . and o 

K min H max 

are given by 


p mi n (h > * W h i> 
p max^ h ^ “ p max^i^ 


exp 


exp 


( 

( 



mm 



H 


max 


h 


h 


) 

) 


(5.2-12) 


where (fu < h < h i+1 ) and the respective scale heights, 
H min and H max' are 9 iven by 


H . 


mm 


h. - h. . 
1 n-1 


in [ p . (h. +l)/p . (h. ) 1 

K mm v i ' K mm v i' J 


(5.2-13) 


H 


max 


h. - h. _ 
i l + l 


in[ Pmax < h i +1 > / Pmax< h i>3 


(5.2-14) 
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A good approximation (neglecting polar motion) for the sat 
ellite height, h, is given by 

h = r - R e (5.2-15) 


where R„ is the mean radius of the Earth, given as 
E 


r e 



Vi - £ e> 


( 2f E - f E ) COS 2 6 


(5.2-16) 


and r = magnitude of the satellite position vector 
= equatorial radius of the Earth 

f E = Earth's flattening coefficient 

= l-(polar radius)/(equatorial radius) 

6 = declination of the satellite (it is assumed that 6 
equals the geocentric latitude of the subsatellite 
point) 

If the density is assumed to be maximum at the apex of the 
bulge, the cosine variation between the maximum and minimum 
density profiles is 


D.(h) = p • (h) 

K 0 V K mm v ' 


[p (h) - 
LK max v ' 


P • (h)] 

K nnn' 7 


cos 


n ( 


(5.2-17) 


where y is the angle between the satellite position vector 
and the apex of the diurnal bulge. 

The cosine function in Equation (5.2—17) can be determined 
directly as 


cos 


n X _ 


1 + COS V 


n/2 


T . tt -jn/2 


U 


B 


2 r 


(5.2-18) 
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where r = satellite position vector 

Ug = unit vector directed toward the apex of the 
diurnal bulge 

For ADEAS, n has been assigned the value of 6. 

The vector Ug has the following components: 


Ug = COS Sg 

X 

cos(a 

s 

+ X) 

(5.2-19) 

U By = cos Sg 

sin(a 

s 

+ X) 

(5.2-20) 

U„ = sin 

B -7 

6 s 


(5.2-21) 


where Sg = declination of the Sun 

ag » right ascension of the Sun 

X = lag angle between the Sun line and the apex of 
the diurnal bulge = 30 degrees 

5. 2. 1.3 Sun and Moon Perturbations 

The perturbations on the spacecraft’s orbit due to its 
attraction by a point mass Sun and/or Moon are optionally 
included in ADEAS. The acceleration experienced by a sat- 
ellite due to the Sun and Moon, expressed in an inertial 
coordinate system, is 



E 



(5.2-22) 


where r = vector from the center of an inertial coordinate 
system to the satellite 

Pk = product of the universal gravitational constant 
and the mass of the kth point mass 

r kp = vector from the kth point mass to the satellite 
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r^p = magnitude of the vector 
k = 1 for the Sun 
k = 2 for the Moon 

In ADEAS, the motion of the spacecraft is referenced to the 
Earth’s position. The desired form for the acceleration is 
obtained by subtracting the acceleration acting on the ref- 
erence body. 



(5.2-23) 


where r k = vector from the kth point mass to the 
body (the Earth) 

r k = magnitude of the vector r k 
from each side of Equation (4-25) to obtain 


reference 



(5 


2-24) 


For the case in which the reference central body, i.e., the 
center of the coordinate system, is the Earth, and the per- 
turbing point masses are the Sun and the Moon, Equa- 
tion (4-27) becomes 


where 


a 


S 




(5.2-26) 
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(5.2-27) 


The quantities and p M are the gravitational constants of 
the Sun and the Moon, respectively, and r c and r M are the 
position vectors of the Sun and the Moon, referenced to the 
Earth. 

5.2. 1.4 Runoe-Kutta Integration 

The orbital equations of motion can be written compactly in 
terms of the 6-component vector 


as 



(5.2-28) 



dX/dt = f (X, t) 


(5.2-29) 


The fourth-order Runge-Kutta numerical integration gives 
values X n of the vector at a sequence of times t , for n = 
1, 2, ..., given the initial value 


x 0 5 X(t 0 ) 


(5.2-30) 


The solution is advanced from time t' to 

n 


^n+l * fc n + h 


(5.2-31) 
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by 


X n+ i = X n + (h/ 6 )(f 1 + 4f 3 + f 4 ) (5.2-32) 


where 


f l E f < X n' 


fc n> 


(5.2-33) 


f_ = f(X + hf,/4, t_ + h/4) (5.2-34) 

2 n 1 n 

f 3 = f(X + hf 2 /2, t n + h/2) (5.2-35) 

f 4 = f(X + h(f x - 2f 2 + 2f 3 ), t n + h) (5.2-36) 

The integration stepsize h is determined by accuracy consid- 
erations. For ADEAS the stepsize is a user input. 

5. 2. 1.5 Hermite Interpolation 

The spacecraft position and velocity are typically required 
at points other than the points Xq, X^, X 2 , ... given by the 

Runge-Kutta integration. For output at intermediate times, 
a Hermite interpolation scheme is used. 

Assume we have the spacecraft position and velocity r^, v^ 
at the n equally-spaced times t. t 'k.+ 2 / ***' fc k+n 

spacing h = t. 1 - . Then, the n-point Hermite interpola- 

tion for the time 


t 


fc k + l + sh 


(5.2-37) 
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is 


r(t) 


v(t) 


where 


g ( i , s) 
y vr v 


with 

A^s) 


n 

- 2 '9rr (i ' 5) 7 k + i + 9 tv< 1 ' s > ’k.i 1 <5 

i = l 

n 

= X) [g vr (i ' s > 7 k.i + V (i ’ s > *k + i ] < 5 

i=l 


g rr (i, s) = [1 - 2C i (s) D i (s) 2 (5 

g rv (i/ s) = C^s) D^s) 2 h (5 

= (2/h)[A i (s) - 2A i (s) B ± C^s) - B i D-^s)] D 

(5 

v (i, s) = [D i (s) + 2 A-^s) C-^s)] D i (s) (5 



n 

j)" 1 (s - k + l)/(i - k) 

k=l 
k^i, j 


(5 
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(5.2-46) 


C^s) = s - i + 1 


n 

D^(s) E [(s - j + l)/(i - j)] (5.2-47) 

j=l 

For the Runge-Kutta integration a four-point interpolator 
(n = 4) is used. A six-point interpolator (n = 6) is used 
for interpolating data from the ephemeris file, as described 
in Section 5.2.2. 

5.2.2 EPHEMERIS FILE 

The position and velocity for a spacecraft at a specified 
time may be evaluated by accessing a standard GTDS EPHEM 
file (Reference 4) . 

5.3 GEOMAGNETIC FIELD VECTOR COMPUTATION 

Whenever the analyst specifies a spacecraft using a magne- 
tometer, the desired magnetic field order must be specified 
between 2 and 10. The computation of the direction of the 
Earth's magnetic field at any time and position is based on 
the spherical harmonic model of Appendix H (Reference 3). 
Specifically, the field is computed using Equations (H-12) 
to (H-15) of that reference. 

5 . 4 SUN AND MOON EPHEMERIS GENERATION 

The following sections specify the computation of the Sun 
and Moon vectors in an Earth-centered coordinate system. 

The spacecraf t-to-Sun and spacecraf t-to-Moon vectors are 
computed by subtracting the Earth-to-spacecraf t vector from 
the Earth-to-Sun and Earth-to-Moon vectors, respectively. . 

5.4.1 SLP FILE 

SLP ephemeris data files are regularly created using a JPL 
planetary ephemeris tape as a data source. 


0450 


5-15 


The solar and lunar positions are calculated from 


^ (li^) ^ 

i=l 


where X(t) = desired position vector 

Ki = ephemeris polynomial coefficients from SLP 
data file 

t = time in seconds from beginning of ephemeris year 

t m = time in seconds from beginning of ephemeris year 
to the midpoint of the curve-fit interval, from 
SLP data file 

N = degree of polynomial plus one 
for Moon, N = 9 
for Sun, N = 5 

5.4.2 ANALYTIC SUN EPHEMERIS COMPUTATION 

The position of the Sun (X , Y , Z ) in mean-of-date GCI co- 

s s s 

ordinates is determined by evaluating a series expansion in 
the longitude of the Sun. The components of the Sun's posi- 
tion are related to the solar orbital elements by 



X s = r s (cos T s ) 


(5.4-2) 


Y s = r s (cos e sin x g ) 


(5.4-3) 


Z s = r s (sin e sin x g ) 

where c = mean obliquity of the ecliptic 

t = solar longitude 
s 

r ■ Earth-to-Sun distance 
s 
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The mean obliquity of the ecliptic (e) is defined as 


c = 23?4522940 - 0?0130125 T c (5.4-5) 

where T represents the full Julian date of the current 
time as the number of Julian centuries since 1900.00: 

( Julian date of c urrent time) - 241502Q t . 0 ( 5 . 4 - 6 ) 

T c ■ 36525.0 

The solar longitude (T g ) i n degrees is defined as 

x s = g s + is + (360°/ir) es sin (as) (5.4-7) 

where G = longitude of perigee of the solar orbit 
s 

= 281?220844 + 0?0000470684 AD 
as = solar mean anomaly 

= 358?475833 + 0?985600267 AD 

es = eccentricity of the Earth's orbit about the Sun 
= 0.01675104 - 0.11444 x 10 -8 AD 

and AD represents the full Julian date of the current time 
as the number of days since 1900.00: 

AD = 36525.0 T, (5.4-8) 

The Earth-to-Sun distance (r g ) is defined as 

— s = 1 + e cos (as) (5.4-9) 

r s S 

5-17 


0450 


where d s = mean solar geocentric distance 
= 149597871.0 kilometers 

5.4.3 ANALYTIC MOON EPHEMERIS COMPUTATION 

The position of the Moon (X . Y . Z 1 in mean-of-date GCI 

m m m 

coordinates is determined by evaluating series expansions in 
the latitude and longitude of the Moon. The components of 
the Moon's position are related to the lunar orbital ele- 
ments by 


X = r (cos 9 cos t ) 
mm m m 


(5.4-10) 


Y m = r m (cos e m sin T n, cos E “ sin e m sin e > (5.4-11) 


m 


m 


m 


Z m = r m (cos 9 m sin T m sin e + sin e m cos e) (5.4-12) 


m 


m 


m 


where e = mean obliquity of the ecliptic 

T m = lunar ecliptic longitude 

9 m = lunar ecliptic latitude 
m 

r = Earth-to-Moon distance 
m 

The mean obliquity of the ecliptic (e) is defined by Equa- 
tions (5.4-5) and (5.4-6). 

The lunar ecliptic longitude (t ) is defined as 

m 

\n - L m + series ! < D m' F m' 4 m' l s> (5.4-13) 


where 


L m = the mean longitude of the 

Moon, measured in the eclip- 
tic plane from the mean 
equinox of date to the mean 
ascending node of the lunar 
orbit, and then along the 
orbit 
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seriesi (D m , F m/ & m , 


The terms within the series 


= 270°26 ' 02 . 99 + (480960?0 

+ 307 O 52’59 , . , 31)T c - 4.08 

2, s ) = a series expansion in the 
angles D m , F m , 2. m , and i s 
presented in Table 5-1 

expansion are defined as 


D = mean elongation of the Moon from the Sun 
m 

_ 350°44'14.95 + (444960?0 + 307°06 ’ 51 . 18)T c (5.4-14) 

" 2 

- 5.17 T c 

argument of latitude of the Moon 

11°15'03.20 + (483120?0 + 82°01 ' 30 . 54)T c (5.4-15) 

" 2 

- 11.56 T£ 

lunar mean anomaly 

296°06 ' 16" 59 + (477000?0 + 198°50 1 56" 79 ) T c (5.4-16) 

" 2 

+ 33.09 T c 

Q, = solar mean anomaly 
s 

= 358°28 , 33 , .'00 + (35640?0 + 359°02 1 59 " 10 ) T c (5.4-17) 



The lunar ecliptic latitude is (© m ) is defined as 


0 m = series 2 (D m , F m , 


m 


a s> 


(5.4-18) 
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Table 5-1 



-0.000607 

sin 

0 

0.011490 

sin 

0 

-0.000267 

sin 

0 

-0.001996 

sin 

0 

-0.000801 

sin 

0 

-0.003238 

sin 

0 

-0.000118 

sin 

0 

0.000138 

sin 

1 

0.000716 

sin 

1 

0.000192 

sin 

1 

-0.000186 

sin 

1 

-0.022236 

sin 

1 

0.10976 

sin 

1 

0.000931 

sin 

1 

-0.000219 

sin 

1 

-0.000999 

sin 

1 

-0.000532 

sin 

1 

-0.000149 

sin 

2 

-0.001026 

sin 

2 

0.003728 

sin 

2 

0.000175 

sin 

3 
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where series 2 ( D m' F i 


m' 


^m' ^s) 


a series expansion in the 
angles D m , F m , S. m , and l s 

presented in Table 5-2 


The Earth-to-Moon distance (r m ) is defined as 



1 + series. 


<v 


F m' V V 


(5.4-19) 


where ^m 

series3 (D m , F m , 2- ra , 2-s) 


= mean lunar geocentric dis- 
tance 

= 384399.06 kilometers 

= a series expansion in the 
angles D m , F m , 2. m , a nd &s 

presented in Table 5—3 


As noted previously, the series are summarized in Tables 5-1 
through 5-3. These tables give the coefficients of each 
trigonometric term and the angles that appear in that term. 

The trigonometric terms contain only integer multiples of 
the four angles. The first term in 0^ is 


= 0,089503 sin (F M ) radians 

The formulas for the lunar position include the leading 
terms of Brown's lunar theory. All perturbation terms with 
amplitudes greater than 50 kilometers are included (21 terms 
in the ecliptic longitude, 11 terms in the ecliptic lati- 
tude, and 11 terms in the distance); this achieves an over- 
all positional accuracy of 1 arc-minute (0.005 radian or 
200 kilometers). 
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Table 5-2. Series for 


Coefficient 
( radians} 


K 

M 

Argument Multiple of 

F D 2. 

M U M 56 

0.089503 

sin 

0 

1 

0 

0 

0.000569 

sin 

0 

1 

2 

0 

-0.003023 

sin 

0 

1 

-2 

0 

-0.000144 

sin 

0 

1 

-2 

1 

0.004897 

sin 

1 

1 

0 

0 

-0.000807 

sin 

1 

1 

-2 

0 

0.004847 

sin 

1 

-1 

0 

0 

-0.000967 

sin 

1 

-1 

-2 

0 

0.000301 

sin 

2 

1 

0 

0 

0.000154 

sin 

2 

-1 

0 

0 

0.000161 

sin 

1 

-1 

2 

0 
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Table 5-3 


Series for 


Coefficient 
( radi a n <3 ^ 


M 

Argument 

f m 

Multiple of 
°M 

a s 

V 1 □ LI 1 Oil J / 

0.0082488 

COS 

0 

0 

2 

0 

0.0005604 

cos 

0 

0 

-2 

1 

0.0003369 

cos 

1 

0 

0 

-1 

-0.0002086 

cos 

1 

-2 

0 

0 

0. 0100247 

cos 

1 

0 

-2 

0 

0.0545008 

cos 

1 

0 

0 

0 

0.0009017 

cos 

1 

0 

2 

0 

0.0004219 

cos 

1 

0 

-2 

1 

-0.0002773 

cos 

1 

0 

0 

1 

0.0029700 

cos 

2 

0 

0 

0 

0.0001817 

cos 

3 

0 

0 

0 
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5.5 STAR CATALOG GENERATION 


ADEAS will enable the user to provide a star catalog through 
user input, by specifying a SKYMAP Run Star Catalog file, or 
by specifying the internal generation of a uniform star dis- 
tribution. 

5.5.1 SKYMAP FILE 

ADEAS will have the capability to access a user-specified 
SKYMAP Run Star Catalog data set (Reference 5) to retrieve 
star magnitudes and positions (x, y, and z components) di- 
rectly from the data set. The information will be contained 
within the data set in a format as specified for the GRO Run 
Star Catalog defined as 


Word 

1 

INTEGER *4 

SKYMAP ID number 

Word 

2 

REAL *4 

x-component of GCI vector, epoch 
2000.0, proper motion to 1990.0 

Word 

3 

REAL *4 

y-component of GCI vector, epoch 
2000.0, proper motion to 1990.0 

Word 

4 

REAL *4 

z-component of GCI vector, epoch 
2000.0, proper motion to 1990.0 

Word 

5 ' 

REAL * 4 

Visual magnitude 

Word 

6 

REAL *4 

Not required 

Word 

7 

REAL * 4 

Not required 

Word 

8 

REAL *4 

Not required 


5.5.2 UNIFORM STAR DISTRIBUTION 

ADEAS will have the capability to distribute n stars uni- 
formly over the celestial sphere, where n is a user input. 
The following method is used to distribute the stars. 

The unit sphere has an area of 4 ir, so if a square patch 
were allocated to each star, it would have a side of 


d = -Ja tr/n radians 


(5.5-1) 
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It is impossible to cover a sphere with nonoverlapping square 
patches, but this equation should be a good approximation to 
the star spacing for large n. 

The uniform distribution algorithm puts the first star at 
the north pole, the second star at the south pole, and any 
remaining stars on m rings of constant latitude, 

where 


m = TRUNC /* (n - 1) (5.5-2) 

TRUNC denotes truncation, i.e., TRUNC(X) is the largest in- 
teger less than or equal to X. This function is chosen to 
give a good distribution for small n and to tend to spacing 
d for large n. The rings are at latitude, 

- (if/ 2 ) [ 1 - (k + l)/(m + 1)] for k = 1, 3, 5, ... (5.5-3) 


\ k = -\_ 1 for k = 2, 4, 6, ... (5.5-4) 


where k runs up to m. This indexing is chosen so that the 
rings move toward the equator for increasing k. 

The algorithm is initialized for n > 2 by 


n . . = n - 2 

remaining 


(5.5-5) 


a 


remaining 


m 


cos X, 


k= 1 


(5.5-6) 
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The Quantity ot • • is the total length around the rings 

xiic remaining 

divided by 2ir. Then, the following computations are per- 
formed for k = 1 to m: 


Set n k = ROUND 


cos X, 


n 




remaining <* remaining / 


(5.5-7) 


where ROUND(X) is the closest integer to X. 

Put n^ stars on the ring at latitude X^ and longitude 2irj/n k 
for j = 0, 1, n^ - 1. Set 


n remaining n remaining n k 


and 


a . = a . — cos X, 

remaining remaining k 

5.6 LINE— OF- SIGHT SENSOR VISIBILITY 

Occultation and interference of all line— of —sight sensors 
will be evaluated as requested by the user. Both computa- 
tions may be independently enabled or disabled for each cen- 
sor. When the occultation and interference computations 
have been enabled for a sensor, the interference computation 
will be evaluated first. 

5.6.1 OCCULTATION 

The occultation of a sensor will be defined as the intersec- 
tion of the line-of-sight vector from the spacecraft to the 
observed object and the disk of the occulting object as de- 
fined by that object's angular radius. Line-of-sight sensors 
may be occulted by the Earth or Moon. 
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Given the position vector of the occulting object (o B ) in 
body coordinates and the radius of this object (r), the an- 
gular radius of this object (9) is computed as: 


0 = 



The sensor is occulted by this object when 


(5.6-1) 


9 > cos 

where Lg = line-of-sight unit vector from the spacecraft to 
the observed object in body coordinates 

5.6.2 INTERFERENCE 

The interference of a sensor, excluding the V-slit and gimbaled 
sensors, will be defined as the intersection of the sensor 
boresight and the disk of the interfering object as defined 
by the sum of that object's angular radius and an interfer- 
ence angle associated with that object as defined by the 
user. The V-slit and gimbaled sensor interference computa- 
tions replace the sensor boresight with the line-of-sight 
vector from the spacecraft to the observed star. Line-of- 
sight sensors may be interfered with by the Earth, Sun, or 
Moon . 

Given the position vector of the interfering object (o 0 ) i- n 
body coordinates and the radius of this object (r), the angu- 
lar radius of this object (0) is computed as defined in Equa- 
tion (5.6-1). The sensor, excluding the V-slit and gimbaled 
sensors, is interfered with by this object when 

(9 + t) > cos 


B 


B 


|o 


B' J 


(5.6-3) 


B 


B 


'B 


(5.6-2) 
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where B = unit boresight vector in body coordinates 

t = user defined interference angle for the interfering 
object being evaluated 

The V-slit and gimbaled sensors, similarly, are interfered 
with by this object when 


(9 + x) > cos 


-1 




B 1 


(5.6-4) 


where Lg = line-of-sight unit vector from the spacecraft to 
the observed object in body coordinates 
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